Addressing closures

[yakir12]

Related to #53, maybe an additional method for optimize that accepts DataType types will help with the closures-speed-bump?
We could then use FastAnonymous and custom made types.

[mlubin]

Do you have a particular example when the overhead of using closures is significant?

[yakir12]

In my case, I'm remapping RGB values using a relationship between semi-random noisy colormaps and some quanta (angle of polarization, degree of polarization, etc.). It means finding what, say, degree of polarization would result in an RGB color that best matches that of a certain pixel. I have 500x500 pixels per image, and a few tera images. I wouldn't bother if the improvement is less than 10 orders of magnitude though.

[johnmyleswhite]

Have you tried making it clear to Julia that the closed over variables are constant when defining your closure? In my experience, that's usually sufficient to remove some of the worst performance costs of closures.

[yakir12]

Here is the version of the code I assume would be the fastest:

function main(a)
    function fun(x)
        # some code that involves a and x
    end
    optimize(fun,-π/2,π/2).minimum
end

Where would I set any constants?
Or in a more closure-like syntax:

function fun(x,a)
    # some code
end
function main(a)
    optimize(x -> fun(x,a),-π/2,π/2).minimum
end

All these examples are very minimal, but I want to understand what you meant with the constants.

[johnmyleswhite]

I believe something like this should get you the perf that you need:

function main(y)
    const a = y
    function fun(x)
        # some code that involves a and x
    end
    optimize(fun,-π/2,π/2).minimum
end

It's worth noting that closures aren't anonymous functions. Your first example is a named function that's also a closure. Your second example is an anonymous function that's also a closure. Being anonymous doesn't make your code more "closure-like".

[yakir12]

Awesome, thanks! I'll confess that I thought that closure and anonymous are the same :P

Didn't do much though, here's a concrete example, RGB is a Color.jl type which is immutable (so const doesn't change much, right?):

using  Optim, Color
function fAoP(b::RGB)
    const a = b
    function f(x)
        x = 6.*x/π + 3. # 0 <= x <= 6
        x < 1 ? 
            1. - a.r + abs(a.g - x) + a.b :
        x < 2 ?
            abs(2. - x - a.r) + 1. - a.g + a.b :
        x < 3 ?
            a.r + 1. - a.g + abs(x - 2. - a.b) :
        x < 4 ?
            a.r + abs(4. - x - a.g) + 1. - a.b :
        x < 5 ?
            abs(x - 4. - a.r)   + a.g + 1. - a.b :
        # x < 6
            1. - a.r + a.g + abs(6. - x -a.b)
    end
    optimize(f,-π/2,π/2).minimum
end

[timholy]

I have 500x500 pixels per image, and a few tera images

Given this use case, you would be distinctly unhappy with FastAnonymous: see the last sentence in https://github.com/timholy/FastAnonymous.jl#inner-workings. You'll create one new type per pixel, and those types are never garbage-collected. And the type-creation itself is going to be very slow, erasing any advantage of using FastAnonymous.

You should be able to use call-overloading to solve this, if Optim can be modified to allow Type inputs in places where Functions are currently required:

julia> using Color

julia> import Base.call

julia> type Yakir{T}
           a::RGB{T}
       end

julia> function call(y::Yakir, x)
           x = 6.*x/π + 3. # 0 <= x <= 6
           x < 1 ? 
               1. - y.a.r + abs(y.a.g - x) + y.a.b :
           x < 2 ?
               abs(2. - x - y.a.r) + 1. - y.a.g + y.a.b :
           x < 3 ?
               y.a.r + 1. - y.a.g + abs(x - 2. - y.a.b) :
           x < 4 ?
               y.a.r + abs(4. - x - y.a.g) + 1. - y.a.b :
           x < 5 ?
               abs(x - 4. - y.a.r)   + y.a.g + 1. - y.a.b :
           # x < 6
               1. - y.a.r + y.a.g + abs(6. - x -y.a.b)
       end
call (generic function with 1004 methods)

julia> b = RGB(rand(), rand(), rand())
Color.RGB{Float64}(0.28964178879399727,0.6508498546848669,0.5720379528319761)

julia> y = Yakir(b)
Yakir{Float64}(Color.RGB{Float64}(0.28964178879399727,0.6508498546848669,0.5720379528319761))

julia> y(3.2)
5.044795833451628

julia> using Optim

WARNING: deprecated syntax "{}" at /home/tim/.julia/v0.4/Calculus/src/symbolic.jl:108.
Use "[]" instead.

WARNING: deprecated syntax "{}" at /home/tim/.julia/v0.4/Calculus/src/symbolic.jl:121.
Use "[]" instead.

WARNING: deprecated syntax "{a=>b, ...}" at /home/tim/.julia/v0.4/Calculus/src/deparse.jl:1.
Use "Dict{Any,Any}(a=>b, ...)" instead.

WARNING: deprecated syntax "{a=>b, ...}" at /home/tim/.julia/v0.4/Optim/src/levenberg_marquardt.jl:41.
Use "Dict{Any,Any}(a=>b, ...)" instead.

WARNING: deprecated syntax "{a=>b, ...}" at /home/tim/.julia/v0.4/Optim/src/levenberg_marquardt.jl:90.
Use "Dict{Any,Any}(a=>b, ...)" instead.

julia> function fAoP(b::RGB)
           f = Yakir(b)
           optimize(f,-π/2,π/2).minimum
       end
fAoP (generic function with 1 method)

julia> fAoP(b)
ERROR: MethodError: `optimize` has no method matching optimize(::Yakir{Float64}, ::Float64, ::Float64)
Closest candidates are:
  optimize{T<:Real}(::Function, ::T<:Real, ::T<:Real)

 in fAoP at none:3

For your use case, the following definition is especially attractive:

function fAoP{T<:RGB}(B::AbstractArray{T})
    f = Yakir(B[1])
    for b in B
        f.a = b
        optimize((f,-π/2,π/2).minimum
    end
end

(which is why I didn't make Yakir immutable).

[yakir12]

Wow, thanks a lot Tim!!! I'll eagerly review this more carefully tomorrow. Again, thank you very much for your time and effort!

[timholy]

In case the point is not apparent:

julia> function fAoP(x, b::RGB)
           const a = b
           function f(x)
               x = 6.*x/π + 3. # 0 <= x <= 6
               x < 1 ? 
                   1. - a.r + abs(a.g - x) + a.b :
               x < 2 ?
                   abs(2. - x - a.r) + 1. - a.g + a.b :
               x < 3 ?
                   a.r + 1. - a.g + abs(x - 2. - a.b) :
               x < 4 ?
                   a.r + abs(4. - x - a.g) + 1. - a.b :
               x < 5 ?
                   abs(x - 4. - a.r)   + a.g + 1. - a.b :
               # x < 6
                   1. - a.r + a.g + abs(6. - x -a.b)
           end
           f(x)
       end
fAoP (generic function with 2 methods)

julia> function callfAoP(B::AbstractArray)
           s = 0.0
           for b in B
               s += fAoP(3.2, b)
           end
           s
       end
callfAoP (generic function with 1 method)

julia> function callcall(B::AbstractArray)
           s = 0.0
           f = Yakir(B[1])
           for b in B
               f.a = b
               s += f(3.2)
           end
           s
       end

julia> B = reinterpret(RGB{Float64}, rand(3, 500, 500), (500, 500));

# After warmup
julia> @time callcall(B)
elapsed time: 0.006607975 seconds (128 bytes allocated)
1.1522790494620719e6

julia> @time callfAoP(B)
elapsed time: 0.203304003 seconds (160 MB allocated, 21.30% gc time in 8 pauses with 0 full sweep)
1.1522790494620719e6

I think many folks will appreciate a 30x speed increase.

[yakir12]

This is awesome! Even if Optim.jl (and for that matter other modules) doesn't implement a Type argument method, I've learned a lot. Really cool.

[timholy]

Certainly Optim can implement "duck typing" for its function inputs. It just needs someone to do the work. Care to give it a whirl?

[yakir12]

"The work" being removing the ::Function from the function call..? This would need to be done for the optimize, brent, and golden_section functions. Very easy.
(I have no idea though if you fork, push request, or what not)

[johnmyleswhite]

Any reason to go to full duck typing rather than use ::Callable annotations?

[timholy]

@johnmyleswhite:

Any reason to go to full duck typing rather than use ::Callable annotations?

Yes:

julia> isa(y, Base.Callable)
false

Callable includes types, but not type-instances. Callable needs to be generalized or eliminated.

@yakir12:

"The work" being removing the ::Function from the function call..?

That, and presumably you need to make things like DifferentiableFunction parametric, i.e.,

immutable DifferentiableFunction{F1,F2,F3}
    f::F1
    g!::F2
    fg!::F3
end

Reason is explained in the FAQ.

[timholy]

Oh, and Yakir: make your changes, run the tests to see if they pass, do a git commit, and then use Pkg.submit("Optim").

[timholy]

Update: FastAnonymous has been rewritten to take advantage of some of the new features of 0.4. Note the main README is for 0.3, but there's a link in the first paragraph for a 0.4-specific README. It now generates real closures, so it has all the features you'd need, @yakir12.

[johnmyleswhite]

Is there a PR for this? If not, I think it's better to close this out in favor of an eventual PR.