# Building an allometric food web - help with matrices

**URL:** <https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709>\
**Category:** New to Julia\
**Created:** [January 16, 2019, 4:40pm UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709 "2019-01-16T16:40:40Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [January 16, 2019, 4:40pm UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709/1 "2019-01-16T16:40:40Z")

</div>

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

```julia
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.

```julia
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:

```julia
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

```julia
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:

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

```

and vector

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

```

Any help is greatly appreciated

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [January 17, 2019, 3:22am UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709/2 "2019-01-17T03:22:05Z")

</div>

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

```julia
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

```julia
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](https://docs.julialang.org/en/v1/manual/arrays/index.html#Broadcasting-1), operating on all elements with a single expression:

```julia
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

```julia
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]`):

```julia
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.

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [January 17, 2019, 12:41pm UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709/3 "2019-01-17T12:41:42Z")

</div>

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

```julia
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

```julia
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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [January 17, 2019, 1:13pm UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709/4 "2019-01-17T13:13:28Z")

</div>

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:

```julia
[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

```

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [January 17, 2019, 4:13pm UTC](https://discourse.julialang.org/t/building-an-allometric-food-web-help-with-matrices/19709/5 "2019-01-17T16:13:43Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/0/0/0022acf92466eb9b15e199fe87586f53f8836a4c.png)

So while it may be possible for me to write

```julia
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?