How to build functions that work efficiently with arguments of different types?

General Usage
, ,
1

For implementing sensitivity analysis, I need to build some functions, for example like the one below.

function just_an_example(a, b, C)
    
    a::Float64
    b::Float64
    C::Vector{Float64}
    return a + b*C
end

(This function does not work at this stage and I am aware of that, I have used for loops to compute the iterations related with the Vector. Btw, as I could recall, while using pythonΒ΄s numpy array, the iterations were automatically solved.)

Is there any efficient way to make the function still usable without modification, when next time other arguments like A and B are Vectors, c is a number instead?
I need you guys open my mind with elegant solutions, thank you!

2 
function just_an_example(a, b, C)
    return a + b*C
end

Julia will compile efficient code for this whenever it is called with different argument types.

1 Like
3

2021-11-15 (2)
2021-11-15 (2)909Γ—678 34.4 KB

I have tried this approach firstly, actually. However it never worked. Could there be any case-specific issue in my problem?

5

You are trying to apply ^ to a vector. You probably want to broadcast it .^. Fortunately, broadcasting also works applied to a scalar.

6

By β€œwork” I think you mean you want it to broadcast. Be aware though that this is not the only meaning possible – Julia takes matrices (and vectors) seriously, and while in this case you get an error that raising a vector to a power doesn’t make sense, you can also silently get answers you may not expect. Raising a matrix to a power for instance does make sense, it means repeated *:

julia> [1 2; 3 4]^2  # == [1 2; 3 4] * [1 2; 3 4]
2Γ—2 Matrix{Int64}:
  7  10
 15  22

julia> [1 2; 3 4].^2
2Γ—2 Matrix{Int64}:
 1   4
 9  16

Often the right solution will be to broadcast the whole function when you call it. But this won’t work so well with un-packing pairs of arguments. You could change your function to avoid that. Or you can include broadcasting inside your function, which will do nothing in the scalar case:

julia> f1(x, (y, z)) = x * y^2 / z;

julia> f1(1, (2, 3))
1.3333333333333333

julia> f1(1, (2, [3,4,5]))  # this already has a meaning
1Γ—3 transpose(::Vector{Float64}) with eltype Float64:
 0.24  0.32  0.4

julia> f1(1, ([2,3,4], 5))  # as in question
ERROR: MethodError: no method matching ^(::Vector{Int64}, ::Int64)

julia> f1.(1, ([2,3,4], 5))  # can't broadcast the whole function
ERROR: BoundsError: attempt to access Int64 at index [2]

julia> f2(x, (y, z)) = @. x * y^2 / z;

julia> f2(1, (2, 3))  # same as f1
1.3333333333333333

julia> f2(1, (2, [3,4,5]))  # NB, this is elementwise ./ not linearsolve /
3-element Vector{Float64}:
 1.3333333333333333
 1.0
 0.8

julia> f2(1, ([2,3,4], 5))
3-element Vector{Float64}:
 0.8
 1.8
 3.2

julia> @code_llvm debuginfo=:none f1(1, (2, 3))  # compare this to next...
define double @julia_f1_9765(i64 signext %0, [2 x i64]* nocapture nonnull readonly align 8 dereferenceable(16) %1) #0 {
top:
  %2 = getelementptr inbounds [2 x i64], [2 x i64]* %1, i64 0, i64 0
  %3 = getelementptr inbounds [2 x i64], [2 x i64]* %1, i64 0, i64 1
  %4 = load i64, i64* %2, align 8
  %5 = mul i64 %4, %0
  %6 = mul i64 %5, %4
  %7 = sitofp i64 %6 to double
  %8 = load i64, i64* %3, align 8
  %9 = sitofp i64 %8 to double
  %10 = fdiv double %7, %9
  ret double %10
}

julia> @code_llvm debuginfo=:none f2(1, (2, 3))  # ... check that broadcast compiles away
define double @julia_f2_9769(i64 signext %0, [2 x i64]* nocapture nonnull readonly align 8 dereferenceable(16) %1) #0 {
top:
  %2 = getelementptr inbounds [2 x i64], [2 x i64]* %1, i64 0, i64 0
  %3 = getelementptr inbounds [2 x i64], [2 x i64]* %1, i64 0, i64 1
  %4 = load i64, i64* %2, align 8
  %5 = mul i64 %4, %0
  %6 = mul i64 %5, %4
  %7 = sitofp i64 %6 to double
  %8 = load i64, i64* %3, align 8
  %9 = sitofp i64 %8 to double
  %10 = fdiv double %7, %9
  ret double %10
}
5 Likes
7

Thank you, this is helpful. However, does this (@.) mean that I have to write the calculation function in one line? What should I do when I want to write a complicated function with multiple lines and still to be able to make the function work, regardless of the type (either a single number or a vector) of the arguments?

8

Can you give an example of such function code?

9

Sure, like the example previously given:
I can write in one line to make the function Stokes_number() work, both when 𝒅_particle is a single number and a vector.

Input:
Stokes_number(π›ˆ, (𝝆_particle, 𝒅_particle), (𝒗_bubble, 𝒅_bubble)) = @. (𝝆_particle * 𝒗_bubble * 𝒅_particle^2) / (9 * π›ˆ * 𝒅_bubble)

test = Stokes_number(0.01, (2.65, 20), (31.6, 0.77*10^-1))

init_vec = [1, 2, 3, 4, 5]
test2= Stokes_number(0.01, (2.65, init_vec), (31.6, 0.77*10^-1))

Output:
Stokes_number (generic function with 1 method)

4.833477633477633e6

5-element Vector{Float64}:
  12083.694083694081
  48334.776334776325
 108753.24675324673
 193339.1053391053
 302092.35209235206

However, when my functions have to be more complex in the future, IΒ΄d prefer to write them in multiple lines, for example:

function Stokes_number(π›ˆ, (𝝆_particle, 𝒅_particle), (𝒗_bubble, 𝒅_bubble))
   ​
   ​𝙆 = (𝝆_particle * 𝒗_bubble * 𝒅_particle^2) / (9 * π›ˆ * 𝒅_bubble)
   ​
   ​return 𝙆
   ​
end

test3 = Stokes_number(0.01, (2.65, 20), (31.6, 0.77*10^-1))
init_vec = [1, 2, 3, 4, 5]
test4 = Stokes_number(0.01, (2.65, init_vec), (31.6, 0.77*10^-1))

Output:
Stokes_number (generic function with 1 method)

4.833477633477633e6

MethodError: no method matching ^(::Vector{Int64}, ::Int64)

This time it wonΒ΄t work. And as @mcabbott mentioned, by coding vector^2 there is a high possibility that something unexpected will come out. I just wanted to know how should I better build my functions, maybe by removing parentheses and avoiding pair-wise arguments?

10 
function Stokes_number(π›ˆ, (𝝆_particle, 𝒅_particle), (𝒗_bubble, 𝒅_bubble))
    K = @. (𝝆_particle * 𝒗_bubble * 𝒅_particle^2) / (9 * π›ˆ * 𝒅_bubble)
    return K
end

or even

function Stokes_number(π›ˆ, (𝝆_particle, 𝒅_particle), (𝒗_bubble, 𝒅_bubble))
    K = @. (𝝆_particle * 𝒗_bubble * 𝒅_particle^2) / 
        (9 * π›ˆ * 𝒅_bubble)
    return K
end

BTW, I got a lot of syntax: invisible character \u200b near column 5 problems when copy-pasting your code. I love unicode characters, but I think they should be used sparingly, and only when there is no good ASCII equivalent. For example, why use 𝙆, 𝒅, 𝒗 when K, d, v is better? Using unicode characters that look almost exactly like ASCII characters, and even have the same meaning, seems like a recipe for trouble.

4 Likes
11

Why return K?

12

Not sure, but the OP wanted to rephrase the function as a multiline function, as an example, so I didn’t change that.

13

I got a lot of syntax: invisible character \u200b near column 5 problems when copy-pasting your code. … why use 𝙆, 𝒅, 𝒗 when K, d, v is better?

Sorry, I didnΒ΄t realize that unicode characters could cause so much trouble. And surely your suggestion is right, I will avoid using them unnecessarily.

Thank you very much for the codes, now the function behaves like IΒ΄ve expected.

1 Like
14

Sorry, I tried your method with further implementations (more complex scenarios). Your method seems only to be able to handle functions with single calculation. In the case of this function below, which needs multiple lines of calculation (for the sake of readability), adding a @. before E0 will cause the E0 not defined problem.

IΒ΄m wondering how experts would build functions in a way, which regardless of the type of one argument (either itΒ΄s a single value or a vector), functions could still be computed automatically. (Of course when two arguments are vector, then itΒ΄s another story, I guess the syntax IΒ΄ve used here will cause more troubles for vector multiplications.)

# The GSE model, Generalized Sutherland Equation model
function CollEffi_GSE(d_particle, d_bubble, 𝝆_particle, 𝝆_fluid, Re)
    E0 = 3 * d_particle / d_bubble
    𝛃 = 12 * d_bubble * 𝝆_fluid / (d_particle * (𝝆_particle - 𝝆_fluid) * Re)
    sin_sqr_term = 2 * 𝛃 * ((1 + 𝛃^2)^0.5 - 𝛃)
    cos_term = sqrt(1 - sin_sqr_term)
    term1_in_parent = 4 * E0 * cos_term * (log(3/E0) - 1.8) / (3 * 𝛃)
    term2_in_parent = 2 * (2/3 + cos_term^3/3 - cos_term) / (𝛃 * sin_sqr_term) 

    Ec = E0 * sin_sqr_term * (exp(term1_in_parent - term2_in_parent))
    return Ec
end

Dp_Vector = collect(1:1:100)
test1 = CollEffi_GSE(Dp_Vector, 0.77*1000, 2.65, 1, 240)
15

@. broadcasts the following complete single expression. If you put it in front of a variable assignment, it converts it to .=, which means in-place updating, and can only work if the assigned-to variable already exists.

In a function body you must add as many dots or @. as needed, in the right places.

In your current function you should add none, though. Instead call CollEffi_GSE itself with a dot, as CollEffi_GSE.().

You should normally not do that. Try to write your functions to accept scalars, and then broadcast the function itself, unless the function inherently is a vector function.

3 Likes
16

IΒ΄ve tried to put a dot before each operator and Base function, according to the documentation of Julia. That was working as expected, however, is obviously not a convenient approach.

Instead call CollEffi_GSE itself with a dot, as CollEffi_GSE.() .

This works so well and I didnΒ΄t realize that my own functions can be called with a dot as well like Base.exp.() etc… Only appreciations to you, friend. :100:

17

Indeed, every function and operator can be called with a dot! It is one of my favourite Julia features. If you’re interested, the blog post that introduces and perhaps best explains the feature is here: More Dots: Syntactic Loop Fusion in Julia

4 Likes