For if statement help


#1

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

#2

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.


#3

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.


#4

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.


#5

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

#6

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.


#7

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.


#8

This would be slightly faster:

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

#10

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

image

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

#11

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.


#12

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

#create a matrix of links:
#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

#13

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