Bridge the gap, please! (from Matlab to Julia)

Sorry a few notes as well,

you shouldn’t benchmark in global scope, you have to put things in functions for benchmark timings to be accurate. You can read about this more in the performance tips in the docs.

Some of the slow timing in your code is because you are re-defining a function over and over again, which causes re-compilation in global scope but not (I think) inside a function. But I’m not sure and you should explore this more.

Here is a benchmark:

julia> struct B{T}
           par::Vector{T}
       end;

julia> function operate(f, a::B)
           f(a.par)
       end;

julia> function foo(n) # you need to benchmark in a function
           for i in 1:n
               if randn() > 0
                   b = B(randn(2))
                   f = sum
               else
                   b = B(rand(3))
                   f = prod
               end
               operate(f, b)
           end
       end;

julia> n = 1e3;

julia> @btime foo($n);
  55.096 μs (1000 allocations: 100.56 KiB)

But to be honest, I’m not 100% sure this is the correct. The variable f is “changing type” based off of a non-deterministic even (the rand call), since all functions are typed on their own. But on the other hand the operate call creates a function barrier so maybe it works itself out. Hopefully someone more informed can chime in.

2 Likes
# create, store and evaluate random functions with random parameters

function as_tuple(n)
    println("as tuple")
    for i = 1:n
        if randn() > 0
            e = (prod, [randn(2)])
        else
            e = (sum, [randn(3)])
        end
        e[1](e[2])
    end
end

function as_array(n)
    println("as array")
    for i = 1:n
        if randn() > 0
            e = [prod, [randn(2)]]
        else
            e = [sum, [randn(3)]]
        end
        e[1](e[2])
    end
end

struct B{T}
    par::Vector{T}
end;
function operate(f, a::B)
    f(a.par)
end;
function as_function(n)
    println("as function")
    for i in 1:n
        if randn() > 0
            f = prod
            b = B(randn(2))
        else
            f = sum
            b = B(rand(3))
        end
        # operate(f, b) # have to add pack/store/unpack
        fb = (f, b)
        operate(fb[1], fb[2]) 
        # fb[1](fb[2]) # does NOT work, why?
    end
end

function run_tests(n)
    @time as_array(n)
    @time as_tuple(n)
    @time as_function(n)
end

run_tests(1e6)

as array
  0.339309 seconds (4.00 M allocations: 297.555 MiB, 7.22% gc time)
as tuple
  0.295214 seconds (3.00 M allocations: 206.005 MiB, 5.65% gc time)
as function
  0.210576 seconds (1.00 M allocations: 99.173 MiB, 7.51% gc time)

My takeaway:
“Think julian” == think in functions, tell the compiler about stable type layout and relations.
Here, we say that f is a function and b (a.par) is a Vector, to be used as its parameters

Still a way to go to understand why

fb[1](fb[2]) 

throws

ERROR: LoadError: MethodError: no method matching iterate(::B{Float64})
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at range.jl:664
  iterate(::Union{LinRange, StepRangeLen}, ::Int64) at range.jl:664
  iterate(::T) where T<:Union{Base.KeySet{var"#s79", var"#s78"} where {var"#s79", var"#s78"<:Dict}, Base.ValueIterator{var"#s77"} where var"#s77"<:Dict} 
at dict.jl:693

It’s because your B struct is not an iterable. It just wraps a container. prod and sum iterate over their argument, so if you want to use them with your struct, you have to implement the iterate interface for your struct.

2 Likes

Should be

fb[1](fb[2].par)

Another way is to define

Base.sum(b::B) = sum(b.par)

Although the the “operate” approach is perfectly fine and more general.

Some resources you will find helpful, if you are transitioning especially from Matlab.

  1. julia quant-econ
  2. techytok
  3. julia wiki
  4. Translate code while making use of discourse, stackoverflow. At first, this annoying as productivity takes a hit but soon you appreciate the move.
1 Like

Thx!

Now happy with my solution below. Suggestions always welcome!

using Base: Float64

# store and evaluate randomly functions with different number of parameters

mutable struct MyStruct
    a::Float64
    b::Float64
end
s1 = MyStruct(1.0, 2.0)

fun1(x) = sin(x)
fun2(x,y) = cos(x) + y
function fun3(x,y,s1)
    s1.b = s1.a*(cos(x) + y)
end
function meval(cmd)
    len = length(cmd)
    if len == 2
        z = cmd[1](cmd[2])
        return
    end
    if len == 3
        z = cmd[1](cmd[2], cmd[3])
        return
    end
    if len == 4
        z = cmd[1](cmd[2], cmd[3], cmd[4])
        return
    end
end

function test1(n)
    for i = 1:n
        if randn() > 0
            cmd = (fun1, randn())
        else
            cmd = (fun3, randn(), randn(), s1)
        end
        meval(cmd)
    end
end

function run_tests(n)
    @time test1(n)
end

run_tests(1e6)

0.330759 seconds (2.01 M allocations: 39.025 MiB, 2.81% gc time, 5.33% compilation time)

Maybe you could tell us what is that you don’t like about the suggestion of @pdeffebach, which is simpler and more general. (Also you are using s1 as a global variable there, which will hurt performance).

1 Like

You don’t have to import Float64, it’s available by default.

1 Like

This was added by Julia or VS Code

Nothing to dislike, I only failed to apply it to the heterogenous case of fun3.
I someone of you succeeds…

using Base: Float64

# create, store and evaluate random functions with random parameters

mutable struct MyStruct
    a::Float64
    b::Float64
end
s1 = MyStruct(1.0, 2.0)

fun1(x) = sin(x)
fun2(x,y) = cos(x) + y
function fun3(x,y,s1)
    s1.b = s1.a*(cos(x) + y)
end
function meval(cmd)
    len = length(cmd)
    if len == 2
        z = cmd[1](cmd[2])
        return
    end
    if len == 3
        z = cmd[1](cmd[2], cmd[3])
        return
    end
    if len == 4
        z = cmd[1](cmd[2], cmd[3], cmd[4])
        return
    end
end

function test1(n)
    for i = 1:n
        if randn() > 0
            cmd = (fun1, randn())
        else
            cmd = (fun3, randn(), randn(), s1)
        end
        meval(cmd)
    end
end

struct B{T}
    par::Vector{T}
end;
function operate(f, a::B)
    f(a.par)
end
function test2(n)
    for i = 1:n
        if randn() > 0
            #cmd = (fun1, randn())
            fb = (fun1, B(randn))
        else
            #cmd = (fun3, randn(), randn(), s1)
            fb = (fun3,  B(randn(), randn(), s1)) 
        end
        #meval(cmd)
        operate(fb[1], fb[2]) 
    end
end

function run_tests(n)
    @time test1(n)
    @time test2(n)
end

run_tests(1e6)

This could be shortened to

function meval(cmd)
    z = cmd[1](cmd[2]...)
end

though I’d return the result directly here instead of modifying a global variable :slight_smile:

2 Likes

So this is Julias version of varargin, thanks.

Repeatedly I hear sceptical remarks about global variables.
How else can I implement in Julia what other languages call objects with properties?
Certainly a typical mental gap of someone coming from another language like Matlab.

In a Matlab event-based simulation I had functions with local variables, all visible to the nested functions.
I then could could use callbacks to these nested functions to modify the states of the objects.
Quite elegant, the only issue being speed. Thats why I want to recreate it in Julia.

@Sukera: your meval replacement throws an error, did you try?

Ah, sorry, I misread your original code - it should be

function meval(cmd)
   z = cmd[1](cmd[2:end]...)
end

instead.

Running, thank you very much!
Not yet familiar with the syntax. :end and … seems telling the same, twice.

I do not know Matlab, so I cannot help you with the equivalence here. But one structure that smells like what you are describing are the “functors”, and the structure of those is that I’ve shown above, as a callable object:

julia> struct Func{T}
         a::T
         b::T
       end

julia> (f::Func)(x,y) = f.a*x + f.b*y

julia> myfunc = Func(1,2)
Func{Int64}(1, 2)

julia> myfunc(5.,10.)
25.0

You can pass that object as a parameter to other functions:

julia> function external(f)
         x = rand()
         y = rand()
         my_result = f(x,y)
         return my_result
       end
external (generic function with 1 method)

julia> external(myfunc)
0.2535047160043371

The cmd[2:end] part is a slicing operation, taking everything from the second index on until the end of the indexable object. f(x...) takes an iterable object x and splats its content into the call to f, assigning it to each successive variable.

It should be mentioned while the syntax is useful, it should be limited to when you know how many elements x will have (and that number is somewhat small).

@lmiq
The example is not consistent and does not run, sorry.

Sorry, while copying/pasting this line, which instantiates myfunc, got lost. (fixed there)

By the way, the most common way to pass functions with variable number of parameters (not variables) to inner functions, is by the use of closures: https://m3g.github.io/JuliaNotes.jl/stable/closures/

julia> const a = 1
1

julia> const b = 2
2

julia> f(x,a,b) = a*x + b*x^2
f (generic function with 1 method)

julia> function outer(f,x)
         result = f(x)
         return result
       end
outer (generic function with 1 method)

julia> outer( x -> f(x,a,b), 5. )
55.0

If you need a variable number of variables, you can again resort to splatting (but this is much less common, and must be used with some care, because it might introduce performance problems if the number of parameters cannot be inferred at compile time):

julia> g(x,y,a,b) = x*a + y*b
g (generic function with 1 method)

julia> function outer2(f,x...)
         result = f(x...)
         return result
       end
outer2 (generic function with 1 method)

julia> outer2( (x,y) -> g(x,y,a,b), 5., 10.)
25.0

julia> outer2( (x) -> f(x,a,b), 5.)
55.0