Too many dots

New to Julia
1

EDIT: In the original version of this post, I included the exact expression for computing the jac output. This allowed other users to run performance tests on dot vs scalar versions of the function. However, I realized that since the expression is based on mathematical equations that have yet to be published, the expression had to be replaced with (EXPRESSION INVOLVING LOTS OF DOTS). Once the work has been published, I can re-include the exact expression and reference the published article.

Hi,

I have a colleague who has coded in Python before and is asking me for advice on a function he has written in Julia. He says that his Julia code, which utilizes the function, runs faster if the input arguments are not type specified. This I can understand. However, to achieve correct functionality, he has coded dot-operations everywhere. The function is:

using NaNMath

function his_function(Îł_b, cosLAB, E_b, m_b, m_t, m_d, m_r, A, B, D)
    # We have to account for the sign in the quadratic solution, before differentiating
    s = -sign.(cos.(Îł_b))

    jac = s .* (EXPRESSION INVOLVING LOTS OF DOTS)

    return jac
end

The input arguments to the function will always have the types:

his_function(Îł_b::Vector{Float64}, cosLAB::Vector{Float64}, E_b::Vector{Float64}, m_b::Float64, m_t::Float64, m_d::Float64, m_r::Float64, A::Vector{Float64}, B::Vector{Float64}, D::Vector{Float64})

My question is simple. I think all the dots make the code look very Python-esque. How would you re-write the code to make it look better, and also run better? Or maybe I am misunderstanding something, and this is fine?

1 Like
2

Perhaps the blog More Dots: Syntactic Loop Fusion in Julia is useful. As mentioned in there, you can use @. in front of the expression to dottify every function call.

If you don’t like broadcasting you can simply use a for loop.

6 Likes
3

Another approach could be to just define the function for scalar values and then use the dot syntax to call it on vectors (or map or loops). This might affect performance but not in the same catastrophic way as in python.

his_function_scalar.(Îł_b, cosLAB, E_b, m_b, m_t, m_d, m_r, A, B, D)
10 Likes
4

Clearly that are lot of terms repeated many times (for instance (E_b .+ 2 * m_b)). It likely that one can get a much more clear and faster code by splitting these calculations in a loop, like:

for i in eachindex(jac, E_b, m_b)
    eb2mb = E_b[i] + 2 * m_b[i]
    ...
    jac[i] = s[i] * eb2mb ...
end

If you put all arrays used inside the eachindex(...) call you safely add @inbounds to the for loop and, likely, achieve close to optimal performance.

1 Like
5

In fact, the performance will probably get better, because the line s = -sign.(cos.(γ_b)) currently allocates two extra temporary arrays that are not needed (one for s and one because the unary - is not dotted). And there might be other dots that were missed, it’s hard to see by inspection.

In fact, I tried this and it got 2x faster on my machine for 1000-element arrays:

julia> Îł_b, cosLAB, E_b, m_b, m_t, m_d, m_r, A, B, D = [rand(1000) for i = 1:10];

julia> @btime his_function($Îł_b, $cosLAB, $E_b, $m_b, $m_t, $m_d, $m_r, $A, $B, $D);
  58.144 ÎĽs (36 allocations: 94.50 KiB)

julia> @btime his_function_nodots.($Îł_b, $cosLAB, $E_b, $m_b, $m_t, $m_d, $m_r, $A, $B, $D);
  28.615 ÎĽs (3 allocations: 7.88 KiB)

If you have a function that logically operates on scalars in Julia, like this one does, there is no reason in general to vectorize it. Just vectorize it when the function is called (if ever! you might also re-write the call site in terms of scalars).

One of the hardest habits to unlearn, coming from Python or Matlab or similar, is the distorted performance model that those languages teach you (“my code is slow, loops are slow, only vectorized library functions are fast”).

8 Likes
6

Thanks for checking. After posting, I also remember some other advantage.
It is also possible to call the scalar function on range and iterators without allocating the inputs and with map! it is also possible to manage when the output vector is allocated.

The main take away here is probably that scalar functions are good in Julia and there is no need to vectorize the internals of the function.

3 Likes
7

Another advantage to making the function scalar and dotting only at the callsite, is that you can then trivially make it work in-place on a pre-allocated buffer, if you want:

jac .= his_function.(...)
2 Likes
8

That’s a bit unfortunate. The old edit is still visible. Can a moderator remove or hide the old edit so the equation disappears? @stevengj?

2 Likes
9

The compiler should factor these out for you with common-subexpression elimination, at least for standard hardware numeric types.

The code might be clearer if you do it manually, of course.

(But you don’t want to do it manually in the original vectorized code because it will result in additional temporary arrays.)

1 Like
10

I don’t see what you are worried about. The equations are incomprehensible without explanation.

(Even if they weren’t, and were a remarkable advance that other people could understand and want to copy, unless you take years to publish it’s not likely that anyone will “scoop” you, and even more unlikely that they would do so without citing you. Lots of code gets published on open-source repositories long before any associated paper is published; I’m not aware of any case where this was a problem.)

Unless these incomprehensible equations are somehow a trade secret?

1 Like
11

True. My concern was merely due to me wishing to exercise good/correct scientific practice.