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