My matrix is actually a vector?


#1

I am missing something as what I had hoped was a matrix seems to be a vector. Continuing from my topic here

My code is

using Distributions
#Constants
Ropt = 100
γ= 2
k = 8.6173324e-5 #Boltzman eV/K
T0 = 273.15 #K
Temp = 293.15 # example of 20 C degrees
q = 1.2 #Hill exponent

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

#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

#create a matrix of links:
#take the link probabilities of the previous matrix and determine if L[i,j] > p
FoodWeb = rand(Uniform(0,1), size(L)...) #what does ... do? This is impossible to search in the Julia documentation
FoodWeb[L .> FoodWeb] .= 1    # Sets F[i,j] = 0  wherever  L[i,j] < F[i,j]
FoodWeb

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)

where

julia> a
50-element Array{Float64,1}:

instead of being a 50x50 matrix. What am I doing incorrectly?


#2

m is a Vector so why is that strange?


#3

It is strange because a is supposed to be a matrix, which shows the attack rate of species with log mass mi on species of log mass mj. It is a matrix in the paper I am using so I must have made a mistake somewhere.

Here is what the code looks like in R.

# generate the log body mass values for each species (log of mi values) 
logBM = runif(50, 10^2, 10^12) # animal body mass 10^2 - 10^12
# sort values 
logBM = sort(logBM)

#create a matrix of link probabilities:
#animals of a certain size MAY eat animals of a particular size, but it is not certain
create.matrix.L <- function(m){ 
  number.of.species = length(m)
  Ropt <-100
  gamma <-2
  matrix.L <- matrix(0, nrow = number.of.species, ncol = number.of.species)
  for(j in 1:number.of.species)
  {
    for (i in 1:number.of.species)              
    {
      matrix.L[i,j] <- (m[i]/(m[j]*Ropt)*exp(1-m[i]/(m[j]*Ropt)))^gamma #Food web Formula from Schneider paper
      if (matrix.L[i,j]<=0.01) #weak links where Lij <= 0.01 were removed 
      {
      matrix.L[i,j] = 0.0
      }
    }
  }
  
  return(matrix.L)
}

# create the L matrix:
L <- create.matrix.L(logBM)

#create a matrix of links:
#take the link probabilities of the previous matrix and determine if L[i,j] > p
create.FW <- function(L){ 
  number.of.species = dim(L)[1]
  mat.of.FW <- matrix(0, nrow = number.of.species, ncol = number.of.species)
  for (j in 1:number.of.species)
  {
    for (i in 1:number.of.species)
    {
      p = runif(1,0,1) #n, min, max
      if (L[i,j]>p) 
      {
        mat.of.FW[i,j] = 1 
      }
    }
  }
  return(mat.of.FW)
}

# create the food web matrix:
food.web <- create.FW(L)

#Constants
k <- 8.6173324e-5 #Boltzman eV/K
T0 <- 273.15 #K
Temp <- 293.15 # example of 20 C degrees

#recreate parameters from paper

create.param.matrix <- function(food.web, m, I, b, c, E){ 
  #FW (to match food web generated), m (bodymass), I (intercept), b (body scaling species i)
  #c (body mass scaling Pred), E (Activation energy)
  number.of.species = length(m)
  param.matrix <- matrix(0, nrow = number.of.species, ncol = number.of.species)
  for (j in 1:ncol(param.matrix))
  {
    for (i in 1:nrow(param.matrix))
    {
      param.matrix[i,j] <- exp(I)*m[i]^b*m[j]^c*exp((E*(T0 - Temp))/(k*T0*Temp)) * food.web[i,j]
    } 
  }
  return(param.matrix)
}

#attack rate a (m^2/s)
matrix.a = create.param.matrix(food.web, logBM, -13.1, 0.25, -0.8, -0.38) #I, b, c, E

#4

In your case size(L) is (50,50).
... used when calling a function will take the elements of the last input parameter and treat them as individual input parameter values.
rand(Uniform(0,1), size(L)...) is interpreted as rand(Uniform(0,1), 50, 50)

a = @. ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38) will apply function ParameterMatrix to the individual elements of m. Therefore the result will be an array.

Instead, consider the following:

function ParameterMatrix(m, I, b, c, E)
    x = @. exp(I)*m^b*m'^c*exp((E*(T0 - Temp))/(k*T0*Temp))
end

a = ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38)

#5

Thank you for your reply, I understand a bit more now. I have changed the code so that it also multiplies with the FoodWeb matrix as I am only interested in the species where L >= FoodWeb

function ParameterMatrix(m, I, b, c, E, 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

I have another function which is supposed to generate a vector this time, but generates an error. I think the error is telling me that I cannot multiply a vector by a matrix, but am not sure.

The error is

K = ParameterVector(m, 10, 0.28, 0.71)
ERROR: MethodError: no method matching ParameterVector(::Array{Float64,1}, ::Int64, ::Float64, ::Float64)
Closest candidates are:
  ParameterVector(::Any, ::Any, ::Any, ::Any, ::Any) at none:3
Stacktrace:
 [1] top-level scope at none:0

And comes from

function ParameterVector(m, I, b, E, 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)

#6

It looks like your call doesn’t match the function: you are giving 4 parameters:

K = ParameterVector(m, 10, 0.28, 0.71)

but it should be 5:

ParameterVector(m, I, b, E, FoodWeb)

#7

Thank you for your reply,
I thought that if something is in the "global environment’ like FoodWeb, it does not have to be specified within the function, similar to

function ParameterMatrix(m, I, b, c, E, FoodWeb)
    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)

#8

Not an expert, but I don’t think it works that way. What you just showed shouldn’t work but I suspect it still does because you once had a function with only 5 arguments (instead of the 6 it has now) and Julia uses the old one because that one matches to your function call.

To achieve what you are looking for you could define the global FoodWeb as the default for the argument:

function ParameterMatrix(m, I, b, c, E, FoodWeb = FoodWeb)

#9

Thank you, I have placed FoodWeb = FoodWeb into the function argument and the error has disappeared.


#10

You’re welcome!

If happen to be new to Julia and you are at all concerned about the performance of your code I would recommend you take a look at the performance tips.
Specifically it would be important to declare your consts as such or give them as parameters to your functions to avoid any globals.