CODEWITHSUNDEEP
julia
Richard Fitzhugh
Richard Fitzhugh

Reputation: 422

Vectorization of Julia expressions involving both variables and arrays

I have the following code that I'm using as part of an attempt to implement the most accurate equation of state for water in Julia 0.6.

struct parameterizedeos
  Tc::Float64
  ρc::Float64
  R::Float64

  #parameters for ideal gas portion
  n₀::Vector{Float64}
  γ₀::Vector{Float64}
end

h2o_n₀ = [-8.3204464837497, 6.6832105275932, 3.00632, 0.012436, 0.97315, 1.27950,
          0.96956, 0.24873]
h2o_γ₀ = [1.28728967, 3.53734222, 7.74073708, 9.24437796, 27.5075105]

function Σ(expr)
    return sum(eval(@. expr))
end

function ig(eos, δ, τ)
  end_ = Σ(eos.n₀[4:8]*log(1-exp(-eos.γ₀)*τ))
  return log(δ) + eos.n₀[1] + eos.n₀[2]*τ + eos.n₀[3]*log(τ) + end_
end

Tc = 647.096
ρc = 322
R = 0.46151805

eos = parameterizedeos(Tc,ρc,R,h2o_n₀,h2o_γ₀)
δ₁ = 838.025/ρc
τ₁ = Tc/500
print(ig(eos,δ₁,τ₁))

Σ is supposed to be a simple form of the corresponding operator from math, while δ and τ use the nomenclature from the linked reference (dimensionless density and temperature respectively). I get LoadError: DimensionMismatch("Cannot multiply two vectors").

I've played around with all sorts of subexamples in the Julia REPL and they all seem to work just as I'd expect. Σ(log(1-exp(-h2o_γ₀)*τ)) vectorizes and sums the elements as expected. Heck, eval(@. h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ) happily returns a 5-element vector. But calling Σ(h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ)) breaks.

I'm a noob at Julia and at arcane things like macros, so anyone could help me figure out what's going on here, that would be great.

Upvotes: 1

Views: 231

Answers (2)

DNF
DNF

Reputation: 12644

The error you get: DimensionMismatch("Cannot multiply two vectors") happens inside the function ig when eos.n₀[4:8]*log(1-exp(-eos.γ₀)*τ) is evaluated, before the result is supposed to be passed to Σ.

Furthermore: @. is expanded at parse time in global scope, it adds dots to function calls and operators. expr is just a variable, not an operator, not a function call, so therefore @. does nothing. eval also does nothing, so sum(eval(@. expr)) is identical to just plain sum(expr). But since the program fails before that point, this is not your problem.

(Edit: Actually eval does more than nothing. It does some work; the results are identical, but the extra work is completely wasted.)

Solution: Delete the function Σ, you don't need it. Re-write the function ig as follows:

function ig(eos, δ, τ)
  end_ = sum(@. eos.n₀[4:8] * log(1 - exp(-eos.γ₀) * τ))
  return log(δ) + eos.n₀[1] + eos.n₀[2]*τ + eos.n₀[3]*log(τ) + end_
end

Edit:

I've played around with all sorts of subexamples in the Julia REPL and they all seem to work just as I'd expect. Σ(log(1-exp(-h2o_γ₀)*τ)) vectorizes and sums the elements as expected. Heck, eval(@. h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ) happily returns a 5-element vector. But calling Σ(h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ)) breaks.

Everything here works as expected. log and exp vectorize automatically, for now, but that behaviour is deprecated. In version 1.0 calling log or exp (or sin or sqrt, etc.) on a vector will be an error.

eval(@. h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ) vectorizes this expression as expected. (But remove the eval, you're just calling it on a value, which it makes no sense to eval.)

And Σ(h2o_n₀[4:8]*log(1-exp(-h2o_γ₀)*τ)) fails because you're trying to multiply two vectors, something that has never worked (you can take the dot-product (using e.g. '* or dot), or do element-wise multiplication, though.) The fact that you later try to pass the result of this into Σis irrelevant.

Upvotes: 1

stefan bachert
stefan bachert

Reputation: 9608

Try this (using .* and log. and exp. instead)

function ig(eos, δ, τ)
  end_ = Σ(eos.n₀[4:8] .* log.(1-exp.(-eos.γ₀)*τ))
  return log(δ) + eos.n₀[1] + eos.n₀[2]*τ + eos.n₀[3]*log(τ) + end_
end

Upvotes: 1

Related Questions