# For if statement help

I would like to create a matrix (FoodWeb) that is the same size as a previous matrix (L).

This second matrix is first filled with 0.
Then each element in L is compared to a random uniform distribution (between 0 and 1) (defined as p in this case).
If the element of L is >= to p, that element in Foodweb becomes 1.

In the end, FoodWeb should only contain 0 or 1 values.

I think my code should look something like this. `i` and `j` are not defined but its the only way I know how to specify a column/row.

``````using Distributions
const Ropt = 100
const Îł= 2
const NumSpecies = 50
const m = sort(rand(Uniform(10^2, 10^12), NumSpecies))

L = @. (m/(m'*Ropt) * exp(1-m/(m'*Ropt)))^Îł
L[L .<= 0.01] .= 0

FoodWeb = fill(0.0,size(L)...)
for FoodWeb[i,:] in 1:NumSpecies
for Foodweb[:,j] in 1:NumSpecies
p = rand(Uniform(0,1), size(L)...)
if L >= p
FoodWeb = 1
end
end
global Foodweb #is this needed?
end
``````

I imagine in a more efficient form it would look something like this, but I do not know how to specify that p is happening within the FoodWeb.

``````FoodWeb = fill(0.0,size(L)...)
p = rand(Uniform(0,1), size(L)...)
Foodweb[L .>= p] .= 1
FoodWeb
``````

If I understand what you want, the following two blocks of code will all accomplish the same thing:

``````FoodWeb = zeros(size(L)...)
p = rand(Uniform(0, 1), size(L)...)
## the explicit loop way
for i in 1:NumSpecies, j in 1:NumSpecies
if L[i, j] >= p[i, j]
FoodWeb[i, j] = 1
end
end
## or, the vectorized way
FoodWeb = zeros(size(L)...)
p = rand(Uniform(0,1), size(L)...)
Foodweb[L .>= p] .= 1
``````

If you donât need the array `p` later on, you could generate the random uniform numbers one at a time inside the loop:

``````for i in 1:NumSpecies, j in 1:NumSpecies
if L[i, j] >= rand(Uniform(0, 1))
FoodWeb[i, j] = 1
end
end
``````

The loop way is a bit more efficient since it doesnât allocate the temporary boolean array `L .>= p`," but if these arrays are only 50x50 and youâre running this once to set up a simulation, the difference in run time is going to be negligible. Write it whichever way is clearer for you to read.

Edit: dot-notation does loop-fusion magic, so thereâs no temporary array created here.

1 Like

Thank you!

Would you be so kind as to explain the difference between the explicit loop and vectorized methods are? That is, I donât understand what it means to be vectorized. I can only see that it doesnât use `[i,j]` notation.

The explicit loop iterates through the indices of each element of the arrays (i.e., the `i'`s and `j`'s), accesses the value stored at each location, and compares them. The vectorized version does the same thing, but Julia writes the loop for you âunder the hood.â

In Matlab, R, and Python (NumPy), for-loops are generally slow. These languages donât have Juliaâs precise type system, so at each iteration of the loop they have to do a bunch of checking to make sure the types of all the variables are compatible and havenât changed since the last iteration, which takes time. To get around this, they implement common operations on vectors/arrays (like +, -, *, /, >, etc.) in a lower-level language like C. When you write `L > p` in Python, Python checks once that the types of `L` and `p` are compatible, then passes them down to a C function. The C function loops through (much faster, since it is a statically typed and compiled language) and compares each element, then passes a boolean array back up to Python. So the standard advice in Matlab, R, and Python is to avoid loops and write things using these âvectorizedâ operations.

In Julia, because of its rigorous type system and JIT compiler, for-loops can run as fast as C, so this is mostly a moot point. You can use vectorized operations or loops and not suffer a performance hit.

Iâll repeat, though, that especially here, where youâre running the operation once in the global scope, the difference in speed is negligible (0.0002 s on my machine), so you should just write it whichever way makes more sense to you. Personally, Iâd use the vectorized versionâitâs more concise and the meaning is clearer to me.

3 Likes

Two tricks to make the vectorised way slightly neater: `zero(L)` is a matrix of the same type and size, but all zero. And you can put the `rand()` into the broadcast like this:

``````foodweb = zero(L)
foodweb[L .> rand.()] .= 1
``````
1 Like

You can also do

``````julia> M = [Int(rand() > p) for i in 1:size(L, 1), j in 1:size(L, 2)]
3Ă3 Array{Int64,2}:
1  1  0
0  0  0
0  0  1
``````

since `rand() > p` itself has a return value of true or false.

1 Like

But I think p is the random number to be generated for each element of L. In which case thereâs an even more compact version:

``````fw = [Int(â > rand()) for â in L]
``````

Since `Bool <: Number` you could possibly omit the `Int`, depending what you do next.

2 Likes

This would be slightly faster:

``````fw = similar(L)
for i in eachindex(L)
fw[i] = ifelse(L[i] >= rand(), 1, 0)
end
``````
1 Like

Continuing my code, I want to build two functions that are based on the formulae

Does the following formulation work?

``````function Feeding(a, B, q, HT)
x = (a*(B')^q)/(1+ sum((HT[:,j:j]) * (a[:,j:j])*B^q))
end

dBasal = function Basal(r,B,K,Feeding)
for i in 1:NumSpecies #part of the equation corresponding to their growth
dBasal[i] = r[i:i,:]*B'*(1-B'/K[:,j:j])
for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
dBasal[j] = dBasal[j] - Feeding(j,i, a, B, q, HT)*B'
end
end
end
``````

Again this is very hard to say without a minimum working example, including definitions of `a`, `b`, `q`, and `HT`.

Some observations on the `Feeding` function:

• You are doing `B^q`, which looks like it does matrix multiplication, whereas your formula looks like every individual element is exponentiated - which would be `B.^q`
• The `j` in `HT[:, j:j]` seems to be undefined, also your sum in the formula runs over all `j` but here youâre picking only one `j`?

In general it might make sense to preallocate a result matrix for `F` of size `i` times `j` and then loop explicitly over `i` and `j` to fill the matrix - in Julia there is no performance penalty for this approach, i.e. you donât have to try and force things into vector/matrix multiplications, and in many cases it helps clarify the code if you write out what happens on an element-by-element basis - it seems to me like yours might be one of the problems where this is the case.

1 Like

Thank you for your comments. You are correct, I should use `B.^q` - I keep forgetting to use `.` when I want it for each element.
For your second point - I become confused with writing scalar and vectorized objects. I know I am mixing the two methods which perhaps isnât possible.

The whole code is

``````using Distributions
#Constants
const Ropt = 100
const Îł= 2
const k = 8.6173324e-5 #Boltzman eV/K
const T0 = 273.15 #K
const Temp = 293.15 # example of 20 C degrees
const q = 1.2 #Hill exponent
const NumSpecies = 50
const B = fill(430.0,50)

const m = sort(rand(Uniform(10^2, 10^12), NumSpecies))

#L is the success curve of consumers aka feeding efficiency, a probability matrix
L = @. (m/(m'*Ropt) * exp(1-m/(m'*Ropt)))^Îł # The @. macro here broadcasts the following scalar expression across `m` and `m'`
# Weak links (Lij <p) are removed
L[L .<= 0.01] .= 0 #  . lets you do things element by element in an array
L

#take the link probabilities of the previous matrix and determine if L[i,j] > p
FoodWeb = zero(L)
FoodWeb[L .> rand.(Uniform(0, 1))] .= 1 # Sets F[i,j] = 0  wherever  L[i,j] < F[i,j]
FoodWeb

function ParameterMatrix(m, I, b, c, E, FoodWeb = FoodWeb)
#FoodWeb (to match food web generated), m (bodymass), I (intercept), b (body scaling species i)
#c (body mass scaling Pred), E (Activation energy)
x = @. (exp(I)*(m^b)*(m'^c)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb
end

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

#handling time HT (s) (TH in Binzer paper)
HT = ParameterMatrix(m,9.66, -0.45, 0.47, 0.26)

function ParameterVector(m, I, b, E, FoodWeb = FoodWeb)
x = @.(exp(I)*(m^b)*exp(E*((T0 - Temp)/(k*T0*Temp)))).*FoodWeb
end
#Carrying capacity K (g/m2)
K = ParameterVector(m, 10, 0.28, 0.71)
#growth rate r (1/s)
r = ParameterVector(m, -15.68, -0.25, -0.84)
#metabolism x ([1/s])
x = ParameterVector(m, -16.54, -0.31, -0.69)
#Feeding Function F

function Feeding(a, B, q, HT)
x = (a*(B')^q)/(1+ sum((HT[:,j:j]) * (a[:,j:j])*B.^q))
end

dBasal = function Basal(r,B,K,Feeding)
for i in 1:NumSpecies #part of the equation corresponding to their growth
dBasal[i] = r[i:i,:]*B'*(1-B'/K[:,j:j])
for j in 1:NumSpecies # part of the equation corresponding to the sum (losses by predation)
dBasal[j] = dBasal[j] - Feeding(j,i, a, B, q, HT)*B'
end
end
end
``````

I think there is something wrong with your 2nd equation â `j` does not appear on the left. If `i` and `j` mean different things then why does `B` appear with both?

But the first is not complicated:

``````julia> a = rand(2,3); B = rand(3); HT = rand(2,3); q=3;

julia> function feed(a,B,q,HT)
Bq = (B') .^ q
low = 1 ./(1 .+ sum(HT .* a .* Bq; dims=2))
a .* Bq .* low
end
feed (generic function with 1 method)
``````

I pulled out `Bq` as powers can be expensive, and doing it inside the broadcast, it would get done 6 times not 3 (with my dimensions) i.e. n^2 not n.

Edit: one more idea is that you may wish to use `Einsum` to stay close to how things are written on paper:

``````julia> using Einsum

julia> function feed2(a,B,q,HT)
Bq = B .^ q
@einsum inner[i] := HT[i,j] * a[i,j] * Bq[j]
@einsum F[i,j] := a[i,j] * Bq[j] /(1 + inner[i])
end
feed2 (generic function with 1 method)

julia> feed(a,B,q,HT) â feed2(a,B,q,HT)
true
``````