I am trying to find out what the going definition for an AbstractArray is in order to properly choose how to dispatch. As usual, this is best understood by an example. Let’s look at two different ways of writing “the same code”.
for i in number_of_iterations
du = f(t,u)
u += dt*du
t += dt
end
# Make a cache for du
for i in number_of_iterations
f!(t,u,du) # update `du` inplace
u .+= dt.*du
t += dt
end
(Yes, for those who know differential equations, this is the Euler method)
If one is doing standard computing and working on arrays, the second method is clearly better because it gets rid of the temporary arrays. However, since I want to provide both support for a broad array of types and good performance, I provide both versions for each method (and there are some good examples which need the first method).
And that’s where the problem comes in. If I only had one method, I would duck type it and leave it for people to see if it works on their weird type. I am kind of thinking that the way to distinguish between the two is by AbstractArray, but I am wondering if that is specific enough. Specifically, the second one will work on (and be faster) any mutable type which implements broadcasting. Is that what an AbstractArray is? Or is there a slightly more general trait I can be dispatching on?
[@dlfivefifty would you mind chiming in on how a Fun fits in? Are they more array-like and should be using the second version in some form?]
Seems easier to just have this decision be made from a user option instead of trying to use dispatch for it. Maybe an inplace = true kw or something.
why not just a version with or without !? (just curious)
An AbstractArray can be read-only. e.g. a Range like 1:3 is an AbstractArray but is immutable. For the methods an array type should implement, see: http://docs.julialang.org/en/latest/manual/arrays.html#Implementation-1
To get a mutable array of the same shape/eltype as a given A::AbstractArray, you use similar(A), or copy!(similar(A), A) to make a mutable copy. In your case, the usual style would be to define both an out-of-place function foo and an in-place function foo!, e.g.:
function foo!(u)
for ...
u .+= ....
end
return u
end
foo(u::AbstractArray) = foo!(copy!(similar(u), u))
foo(u) = foo!(collect(u)) # handle arbitrary iterable u
3 Likes
That last version, with defining the foo! type and then another version that does foo on a copy seems really ideomatic in julia code, at least it is a design pattern I have encountered numerous times.
Is the above way of coding this functionality the preferred / most ideomatic? I am wondering whether just a short function mcopy(u) = copy!(similar(u, u)) would be useful, or even a macro that could be applied to foo! that would produce the relevant foos? I could definitely see myself tripping over not realizing that an array may be immutable, i.e. I have expressions like foo(u) = foo!(copy(u)) many places in my code.
That really convinces me that I was trying to be too fancy and should choose the style based off of the user input. I already place a type-parameter inplace on the problem definition when it’s created (by counting args of the dispatch with the most args, the very much so most common case is to have only one dispatch or have many which add things like Jacobians so that heuristic works very well here). As @kristoffer.carlsson suggests, I should just use that to dispatch instead of some other criteria. Normally I would use something like @stevengj’s proposal, but in this case I want to guerentee that output types match the user’s input types which that violates (but in most places, the convenience is good).
Also, to truly do this for AbstractArrays, you need something more than copy!. I have recursivecopy! in RecursiveArrayTools.jl for doing something like this, because the above breaks on Vector{Vector} etc. I guess that’s another reason why users should be able to choose the form themselves.
Thanks for steering me in the new direction. I think I got stuck for a bit on one silly idea!
Actually, there is a function Base.copymutable(u) that does exactly this, with the bonus that it gives collect(u) for generic iterable types. (See fix and consolidate uses of copy(a) and copy(similar(a),a) in to copymutable(a) by stevengj · Pull Request #16620 · JuliaLang/julia · GitHub) Maybe it should be exported?
1 Like
Funs should work with either form. They are not “array-like” however as there is no getindex.
(Originally, as in chebfun, they were treated as continuously indexed vectors, so getindex was used for function evaluation, and .* was mandatory for multiplication. This analogy was dropped as it wasn’t very intuitive and overloading array interface routines caused issues.)
If (as suggested I think?) that a preferred design pattern is to generally aim to define an inplace foo!(u) and then define a foo(u) = foo!(copymutable(u)), then yes, it would be a great convenience if that was exported 
Oh okay, so they don’t index but still broadcast? Anything that can index can broadcast, right? So boradcast is strictly more general?
[Yeah, the decision is to choose the form via the user given function. This is too subtle to accuracy predict.]
Correct. The support for broadcasting was thrown in per @stevengj’s request, so probably a bit more thought could go into it. For example, boh the broadcasting and non-broadcasting syntax (exp(f) vs exp.(f)) mean the same thing though they use different algorithms.