Use similar in cache construction

[lstagner]

I have been writing a vector calculus package that uses ForwardDiff

With this package you can do things like

grad(phi::Function, x::Coordinate) #Calculate gradient of scalar field phi
laplacian(phi, x::Coordinate) #Calculate laplacian of of scalar field phi
div(A::Function, x::Coordinate) #Calculate divergence of vector field A
curl(A::Function, x::Coordinate) #Calculate curl (Doesn't work yet)

These routines (should) work in any arbitrary coordinate system, even non-orthogonal curvilinear systems.
This is possible because of the Coordinate type carries around a transformation function.

type Coordinate{T,F<:Function} <: AbstractArray{T,1}
    u::Vector{T} #Coordinates
    R::F         #Transformation (x,y,z,...) = R(u₁,u₂,u₃,...)
end

So when you creates co/contravariant vectors it also calculates the correct metric.

The way the Jacobian/Hessian cache was constructed in ForwardDiff was preventing me from using the Coordinate type. Changing the cache constructors to use similar allows me to use the Coordinate type with some caveats.

I have been unable to resolve an issue when usecache=true is set. It hits the type assertion in multithread_(jacobian/hessian)_cachefetch!. As far as I can tell this shouldn't happen. This also cause the "miscellaneous" tests to fail.

For now I have been setting usecache=false and everything works as expected.

[jrevels]

Hi, thanks for the contribution!

As it turns out, I would love to generalize the cache to arbitrary array types, but the fact that we need a manual type assertion here makes it tough, and we can't get rid of it since it's required to combat the type instability from the Dict access. The use of similar in the cache constructors is an attractive solution - I've wanted to switch it to that for a while. The problem is that we can't use that solution unless the type assertions are done properly. Constructing the array via similar in order to retrieve the desired type technically yields the correct result, but is a bit extreme. I'm hoping that we can solve this problem in a way that doesn't cost an array allocation.

Pulling from my comment here, if it's any help:

It would be nice if there were something like similar_type we could call on array-likes, which would allow us to generalize the array types used by the cache types.

[lstagner]

So we need similar_type.

Perhaps we can adapt what the people over at StaticArrays.jl have done.

@andyferris Do you think it would be difficult to generalize StaticArrays similar_type?

[andyferris]

Should be easy enough. What exactly are you proposing? Adding similar_type as a requirement to your inputs?

I don't see why similar_type and similar should only apply to static arrays or arrays in general, respectively. Make sure the first argument is the thing you want to be similar to, and you won't have any clashes with the StaticArray interface. We might also want similar_type added to Base for general AbstractArrays, but we need to have a much larger discussion about the future of the AbstractArray interface for that.

similar_type has origins in FixedSizeArrays.jl, so I'll ping @SimonDanisch, @c42f and @timholy for completeness.

[lstagner]

We basically need something that will swap out the eltype of an AbstractArray and replace it with a Dual type.

Something like this (but not this)

@generated function swap_eltype{T,S}(x::AbstractArray{T}, ::Type{S})
    tname = x.name.name
    tparams = x.parameters
    newtype = :($(tname){$([t == T ? S : t for t in tparams]...)})
    return newtype
end
swap_eltype(rand(Float64,2,2), Int64) == Array{Int64,2}
true

This works but it has a lot of issues (repeated types, ...) I am not sure if meta-programming is even the right solution considering the all the different combinations of type parameters. It would be like whack-a-mole trying to handle all the different variations.

I think the best solution would be to people define similar_type like function for the abstract array they use.

[andyferris]

There is a "safe" way of replacing the eltype by using myType.name.primary and supertype, implemented in both FixedSizeArrays and StaticArrays. The idea is to trace how T in AbstractArray{T,N} descends through the type tree. If the array type cannot satisfy the new eltype, similar_type could even revert to Array for you.

My preference would be to have this defined in Base, but given that that won't be until v0.6, I guess we might as well put it in some package(s). I could add the general definition to StaticArrays but then you'd have to depend on it... on the other hand if it appears in multiple packages that's a pain. If it's just used internally, then give it a new (unexported) name like swap_eltype and let people overload that when it's strictly necessary (i.e. almost never).

[lstagner]

Thanks for the advice. I was able to get something that should work for any AbstractArray

function eltype_param_number(T)
    if T.name.name == :AbstractArray
        return 1
    else
        T_super = supertype(T)
        param_number = eltype_param_number(T_super)
        tv = T_super.parameters[param_number]
        for i = 1:T.parameters.length
            if tv == T.parameters[i]
                return i
            end
        end
    end
end

@generated function replace_eltype{T,S}(x::AbstractArray{T}, ::Type{S})
    pnum = eltype_param_number(x.name.primary)
    tname = x.name.name
    tparams = collect(x.parameters)
    tparams[pnum] = S
    newtype = :($(tname){$(tparams...)})
    return newtype
end

julia> replace_eltype(rand(Float64,2,2), ForwardDiff.Dual{1,Float64})
Array{ForwardDiff.Dual{1,Float64},2}

@jrevels Would you be OK with this?

[andyferris]

Looks good to me. Ideally this is a @pure function but I did encounter a crash once with that...

[lstagner]

OK I updated the branch to use replace_eltype. It also broke my VectorCalculus package in a weird way.

In my tests I get the following error

Differential Operations: Error During Test
  Got an exception of type UndefVarError outside of a @test
  UndefVarError: Coordinate not defined
   at multithread_jacobian_cachefetch!(::VectorCalculus.Coordinate{Float64,VectorCalculus.#Cartesian}, ::ForwardDiff.Chunk{3}, ::Bool, ::Bool) at cache.jl:60
   at jacobian_cachefetch! at cache.jl:63 [inlined]
...

which is weird since it (Coordinate) is defined.

[lstagner]

The error I see with my package also occurs with other packages as well

julia> using ForwardDiff
julia> using StaticArrays
julia> x = @MVector [1,2,3]
3-element StaticArrays.MVector{3,Int64}:
 1
 2
 3
julia> f(x) = prod(x)
f (generic function with 1 method)
julia> ForwardDiff.gradient(f,x)
------------------------------------------------------------------------------------------
UndefVarError                                           Stacktrace (most recent call last)
[#6] — gradient(::Function, ::StaticArrays.MVector{3,Int64})
       ⌙ at gradient.jl:22 (repeats 2 times)
[#5] — #gradient#33(::Bool, ::Bool, ::Function, ::Function, ::StaticArrays.MVector{3,Int64}, ::ForwardDiff.Chunk{3})
       ⌙ at gradient.jl:23
[#4] — vector_mode_gradient
       ⌙ at gradient.jl:94 [inlined]
[#3] — compute_vector_mode_gradient(::#f, ::StaticArrays.MVector{3,Int64}, ::ForwardDiff.Chunk{3}, ::Bool)
       ⌙ at gradient.jl:87
[#2] — jacobian_cachefetch!
       ⌙ at cache.jl:63 [inlined]
[#1] — multithread_jacobian_cachefetch!(::StaticArrays.MVector{3,Int64}, ::ForwardDiff.Chunk{3}, ::Bool, ::Bool)
       ⌙ at cache.jl:60
UndefVarError: MVector not defined

It's like it forgets that the module is loaded?

@jrevels @andyferris Any idea what could cause this?

Edit: Heres what Gallium gives

1|debug > n
In /home/lstagner/.julia/v0.5/ForwardDiff/src/cache.jl:52
59      end
60      return result::NTuple{$NTHREADS,JacobianCache{N,S,jacobian_dual_type(x, chunk)}}
61  end
62  

About to run: (ForwardDiff.jacobian_dual_type)([1,2,3],ForwardDiff.Chunk{3}())
1|debug > s
In /home/lstagner/.julia/v0.5/ForwardDiff/src/cache.jl:46
45  
46  @inline jacobian_dual_type{T,M,N}(arr::AbstractArray{T,M}, ::Chunk{N}) = replace_eltype(arr, Dual{N,T})
47  
48  Base.copy(cache::JacobianCache) = JacobianCache(copy(cache.duals), cache.seeds)

About to run: (Core.apply_type)(ForwardDiff.Dual{N,T<:Real},3,Int64)
1|julia > arr
3-element StaticArrays.MVector{3,Int64}:
 1
 2
 3


1|debug > s
------------------------------------------------------------------------------------------
UndefVarError                                           Stacktrace (most recent call last)

It seems to be some issue with apply_type

[lstagner]

Ha! Take that stupid bug.

It's all good now. 😄

[jrevels]

Thanks for doing all this work, the type tree traversal stuff is pretty interesting. The main downside is that this PR, as it stands, would be a breaking change, since the type assertion will error out for non-AbstractArray inputs. Ref #107, the issue that prompted us to add support for non-AbstractArray inputs. Could we incorporate this work without sacrificing that support?

[lstagner]

I added a fallback method for replace_eltype for non-AbstractArray types. However it would require that whatever type had a method for similar that returns a mutable array-like object. Unfortunately FixedSizeArrays doesn't have similar defined.

[lstagner]

I should also mention that this also requires that the AbstractArray passed to ForwardDiff is mutable.

Added fallback in latest commit

[jrevels]

Unfortunately FixedSizeArrays doesn't have similar defined.

Yeah, it probably doesn't make much sense to define similar on an immutable array.

I should also mention that this also requires that the AbstractArray passed to ForwardDiff is mutable.

That's probably okay. I just played tried it out and realized that ForwardDiff didn't support FixedSizeVectors anyway (I thought it did, but apparently not). At some point, it would be cool to truly support immutable inputs, at least for vector-mode, but that's out of the scope of this PR.

It's pretty unfortunate that the fallback here requires an allocation - how about making the fallback a @generated function, use similar to construct a zero-element array with the same dimensions, and then return the type? That way, the allocation is at least tiny and only occurs the first time the function is called.

[lstagner]

I don't think that will work because it assumes that the size of the array is not a part of its type.

Using something like

julia> @generated function replace_eltype{T}(x, ::Type{T})
           return :(typeof(similar(x, $(T), (repeated(0,length(x))...))))
       end
replace_eltype (generic function with 1 method)

on a StaticArray gives

julia> using StaticArrays

julia> svec = @SVector rand(3)
3-element StaticArrays.SVector{3,Float64}:
 0.202534 
 0.0077592
 0.401944 

julia> replace_eltype(svec, Int64)
StaticArrays.MArray{(0,0,0),Int64,3,0}

To do this you would have to know whether the type had size parameters in advance.
We wont be able to handle every different type in this way. The only reason we can do it for AbstractArrays is because it has a known type structure.

Until we have something like similar_type in Base I think this is the best we can do.

[andyferris]

As a comment to all of the above, it really bums me that ForwardDiff doesn't support immutable containers. My use case involves small vectors (like points in 3D) and using SVector{3,Dual{Float64,...}} and manually unrwrapping is signficantly faster than waiting for that allocation!

Could we use map and so on to add and remove the dual parts? Then the container type doesn't matter. It would work for SVector{3}, standard Array, ArrayFire AFArray, distributed arrays, and other non-array containers too (including FixedSizeArrays)

[jrevels]

@lstagner Hmm, I didn't think of that. That's unfortunate. With that in mind, it might be best to only support AbstractArrays out-of-the-box, and tell users that their array-like type has to overload ForwardDiff.replace_eltype if they want it to work. I think that would be more admissible than adding in an additional array allocation; there's not much point to the cache if subsequent accesses end up consuming a bunch of memory.

@andyferris Vector-mode support for small immutables would be useful, but like I said, it's outside of the scope of this PR. Feel free to open an issue or PR.

[lstagner]

it might be best to only support AbstractArrays out-of-the-box, and tell users that their array-like type has to overload ForwardDiff.replace_eltype if they want it to work.

Do you want me to remove the replace_eltype fallback then?

[jrevels]

Do you want me to remove the replace_eltype fallback then?

Yeah, I think that would be best. Thanks again for all this work!

[andyferris]

@andyferris Vector-mode support for small immutables would be useful, but like I said, it's outside of the scope of this PR. Feel free to open an issue or PR.

Understood, I may just do that one day :) But this idea isn't just for "small" immutables: map is wonderful as it operates optimally from the smallest SVector to the largest DistributedArray and everything in between.

[lstagner]

Merge when ready.

[jrevels]

Can you change the sparse test to reflect the new behavior?

Also it would be good to run the benchmarks to make sure this doesn't incur any performance regressions.

[lstagner]

I updated the tests and here is the benchmark file.
results_julia_v05_forwarddiff_ls_similarcache.jld.zip

Heres my version info I ran the benchmark with

julia> versioninfo()
Julia Version 0.5.0-rc3+0
Commit e6f843b* (2016-08-22 23:43 UTC)
Platform Info:
  System: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, broadwell)

[lstagner]

Hmm do I need to worry about the v0.4 failure?

[jrevels]

Hmm do I need to worry about the v0.4 failure?

It should be easy to fix. Instead of doing @inline @generated function..., you can splice Expr(:meta, :inline) into the returned expression.

[lstagner]

OK tests are passing now.

[jrevels]

Awesome! I ran the benchmarks here versus master, and didn't encounter any reproducible regressions. I'll tag a release with this change after I update some docs later.

Thanks again for this, @lstagner!