exploit fusing broadcasts in Julia 0.5+

[stevengj]

If you implement broadcast(op, f::Fun) = Fun(x -> f(x), domain(f)) then you will be able to do e.g. sin.(erf.(f)) where x::Fun and it will automatically turn into a single fused Fun(x -> sin(erf(f(x))), domain(f)) constructor call.

Furthermore, you will probably want broadcast{T<:Union{Number,Fun}}(op, f::T...) = ... so that you can do operations combining multiple Fun objects and numbers, e.g. atan2.(f,g) or besselj.(3, f).

(In 0.6, this will automatically give you fusing Fun constructors for binary operators like .+ and .^, too.)

See JuliaLang/julia#16285

[stevengj]

You can also implement broadcast! to get in-place operations on Fun, e.g. f .= sin(g) to over-write the coefficient array of f in-place.

[dlfivefifty]

I guess the idea is to get all functions without having to override them?

This will work for smooth functions. I guess the existing overrides that take into account singularities can remain, without the . Notation.

Sent from my iPad

On 13 Aug 2016, at 5:30 AM, Steven G. Johnson notifications@github.com wrote:

If you implement broadcast(op, f::Fun) = Fun(x -> f(x), domain(f)) then you will be able to do e.g. sin.(erf.(f)) where x::Fun and it will automatically turn into a single fused Fun(x -> sin(erf(f(x))), domain(f)) constructor call.

Furthermore, you will probably want broadcast{T<:Union{Number,Fun}}(op, f::T...) = ... so that you can do operations combining multiple Fun objects and numbers, e.g. atan2.(f,g) or besselj.(3, f).

(In 0.6, this will automatically give you fusing Fun constructors for binary operators like .+ and .^, too.)

See JuliaLang/julia#16285

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub, or mute the thread.

[stevengj]

Yes. Also, you avoid having to construct the Chebyshev coefficients multiple times in nested calls.

[stevengj]

You may also want to override broadcast for Domain etc., so that you can do e.g. f = sin.(Domain([0,1])).

[dlfivefifty]

I guess that makes sense if a domain is identified by the identity function on the domain, which is consistent with the Fun(Domain([-1,1])) notation.

Though I'm thinking of changing [-1,1] to -1..1 to mean the unit interval.

[dlfivefifty]

Is there a natural way to do singularity detection a la automatic differentiation? We can determine the range of a function automatically, so can represent it as an Interval, so for example if the range of f is 0..1 and then sqrt.(f) is called, perhaps the sqrt decay can be detected automatically

[stevengj]

Can't you detect the presence (if not the location) of singularities by noticing that the coefficients are not decaying exponentially?

[stevengj]

Note that .. is parsed as an operator, so you can define ..(a,b) = Domain([a,b]) and then e.g. 3..4 defines a domain.

[dlfivefifty]

I don't think looking at the coefficients will be robust enough.

Sent from my iPhone

On 13 Aug 2016, at 12:00, Steven G. Johnson notifications@github.com wrote:

Can't you detect the presence (if not the location) of singularities by noticing that the coefficients are not decaying exponentially?

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub, or mute the thread.

[dlfivefifty]

A basic version of this is now implemented in the development branch.

[dlfivefifty]

@stevengj If there is a vector-valued Fun, how would you expect broadcasting to work? For example, which line should be true:

f = Fun(x ->[exp(x),sin(x)])
norm.(f) == Fun(x ->norm([exp(x),sin(x)]))  # broadcast over x
norm.(f) == Fun(x ->[norm(exp(x)),norm(sin(x))])  # broadcast over vector index

[stevengj]

I would expect norm.(f) == Fun(x ->norm([exp(x),sin(x)])). Think of a Fun as an "array" of values indexed by x.

[dlfivefifty]

I agree, and that is what has been implemented in "development": norm.(::Fun) broadcasts over the domain. Note that arrays of Funs and array-valued Funs behave differently:

f = Fun(x ->[exp(x),sin(x)])
F = [Fun(exp),Fun(sin)]
norm.(f) == Fun(x ->norm([exp(x),sin(x)]))
norm.(F) == [norm(Fun(exp)),norm(Fun(sin))] # isa Vector{Float64}

When you simultaneously broadcast with both Fun and Array arguments, the Array takes precedence:

norm.(f,[1,2]) == norm.(F,[1,2]) == [norm(f[1],1),norm(f[2],2)] # isa Vector{Float64}

Converting the Array to a Fun broadcasts over x again:

norm.(f,Fun([1,2])) == Fun(x->[norm(f[1](x),1),norm(f[2](x),2)]) # isa Fun

EDIT: Ignore the last one.