Hi,
I want to create a specialised high-performance automatic differentiation type. So I propose to use StaticArrays, so that temporary variables in mathematical expressions do not result in pointer allocations, deallocation, and garbage collection. It looks like that:
immutable Adiff{nd}<: Real
x ::Float64
dx ::SVector{nd,Float64}
function Adiff(x::Float64,n::Integer,i::Integer)
this = new(x, allocate space for n Float64, pretty please )
this.dx[:] = 0.
this.dx[i] = 1.
end
end
You will note how I try to charm Julia. Doesn’t help. How do I tell the language that the user said “n”, so I want dx to be an array of size n. And of course, this has to happen in such a way that, when I create an array of Adiffs (in a factory function), the compiler knows its an array of identicaly sized elements.
The way you’ve written this won’t quite work (or at least, it won’t be fast). The problem is that the type of your Adiff depends on nd, and you’re asking that to depend on the value of n in the Adiff(x, n, i) constructor. That means that Julia won’t be able to infer whether Adiff(x, n, i) returns an Adiff{1}, an Adiff{2}, etc. So that will be type-unstable and slow.
Fortunately, there’s a better way. When you call the Adiff constructor, just use the {} syntax to provide nd directly:
immutable Adiff{N} <: Real
x::Float64
dx::SVector{N, Float64}
function Adiff{N}(x::Float64, i::Integer) where {N}
new(x, setindex(zeros(SVector{N, Float64}), 1, i))
end
end
julia> Adiff{2}(1.0, 2)
Adiff{2}(1.0, [0.0, 1.0])
julia> Adiff{2}(1.0, 1)
Adiff{2}(1.0, [1.0, 0.0])
julia> Adiff{5}(1.0, 4)
Adiff{5}(1.0, [0.0, 0.0, 0.0, 1.0, 0.0])
Note the use of setindex to take an all-zeros SVector and return a new SVector with one element set to 1. That’s because SVectors are immutable, so you can’t do dx[i] = 1.
Just so you’re aware, ForwardDiff exists, and is awesome and fairly mature at this point. You could try to use that first before writing your own. Even if you don’t, you might at least get some inspiration from what the author of that package did. See e.g. ForwardDiff.Dual, which is analogous to the ADiff type you’re trying to write.
Excellent. I have a lot to learn to understand this precious snippet of code Adiff{N} OK a parameterised constructor for a parameterised type, makes sense. where {N} where do I find that in the manual? setindex(zeros(SVector{N, Float64}), 1, i) I understand that my type being immutable, this.dx[i] = 1. won’t work. But then if SVector is immutable too, neither should setindex unless setindex does not work in place.
I try to run your code and get “WARNING: static parameter N does not occur in signature for Type at C:\Users\philippem\home\GIT\julia\modules\adiff.jl:10” I am using Julia 0.5 (JuliaPro), is that the problem?
And that takes me to another unformulated question: can I construct my type without involving a copying operation? My impression is that neither the code you or Tamas_Papp works in place. Should I worry? Could I use a generator (i==j ? 1. : 0. for j == 1:N) ?
Thanks tkoolen. Yes, I am aware. I find that implementing dual numbers is a great way to get into a programming language. It never fails to raise a lot of interesting questions. : )
Ah, sorry, yes, my snipped was v0.6 syntax. Happy to help if you want to port this to v0.5, but I’d suggest moving to v0.6 as soon as possible anyway, since it’s got some nice improvements and it’s much closer to what Julia 1.0 will look like.
Re: setindex() note that this is actually a different function than the regular setindex!() (see the !). setindex!() is the in-place mutating method that is called when you do x[i] = 1.
setindex() is a function which takes the same arguments but returns a new StaticArray with the given element set, so it doesn’t operate in place. You’re right to be concerned about creating copies in general, but in this case it’s fine. StaticArrays types are, in general, going to be stored on the stack very efficiently, and LLVM is very good at optimizing things like this setindex() trick. So, at worst, setindex() will do a very cheap copy on the stack (and thus not allocate any memory on the heap), but it’s likely that LLVM can do even better and eliminate the copy. So feel free to setindex() at will.
)