Partial derivatives of TaylorN (#164)

* Add arbitrary partial derivatives of TaylorN

* Add tests and docs

* Add a new method of derivative, and use same convention as Taylor1's derivative

* Add methods for getcoeff involving tuples

docs/src/examples.md

@@ -112,7 +112,7 @@ We note that the above functions use expansions in `Int128`. This is actually
 required, since some coefficients are larger than `typemax(Int)`:
 
 ```@repl fateman
-getcoeff(f2, [1,6,7,20]) # coefficient of x y^6 z^7 w^{20}
+getcoeff(f2, (1,6,7,20)) # coefficient of x y^6 z^7 w^{20}
 ans > typemax(Int)
 length(f2)
 sum(TaylorSeries.size_table)

docs/src/userguide.md

@@ -304,11 +304,12 @@ exy = exp(x+y)
 
 The function [`getcoeff`](@ref)
 gives the normalized coefficient of the polynomial that corresponds to the
-monomial specified by a vector `v` containing the powers. For instance, for
+monomial specified by the tuple or vector `v` containing the powers.
+For instance, for
 the polynomial `exy` above, the coefficient of the monomial ``x^3 y^5`` is
-
+obtained using `getcoeff(exy, (3,5))` or `getcoeff(exy, [3,5])`.
 ```@repl userguide
-getcoeff(exy, [3,5])
+getcoeff(exy, (3,5))
 rationalize(ans)
 ```
 
@@ -344,6 +345,21 @@ an error is thrown.
 derivative( q, 3 )   # error, since we are dealing with 2 variables
 ```
 
+To obtain more specific partial derivatives we have two specialized methods
+that involve a tuple, which represents the number of derivatives with
+respect to each variable (so the tuple's length has to be the
+same as the actual number of variables). These methods either return
+the `TaylorN` object in question, or the coefficient corresponding to
+the specified tuple, normalized by the factorials defined by the tuple.
+The latter is in essence the 0-th order coefficient of the former.
+
+```@repl userguide
+derivative(p, (2,1)) # two derivatives on :x and one on :y
+derivative((2,1), p) # 0-th order coefficient of the previous expression
+derivative(p, (1,1)) # one derivative on :x and one on :y
+derivative((1,1), p) # 0-th order coefficient of the previous expression
+```
+
 Integration with respect to the `r`-th variable for
 `HomogeneousPolynomial`s and `TaylorN` objects is obtained
 using [`integrate`](@ref). Note that `integrate` for `TaylorN`

src/auxiliary.jl

@@ -94,15 +94,17 @@ setindex!(a::Taylor1{T}, x::Array{T,1}, c::Colon) where {T<:Number} = a.coeffs[c
 """
     getcoeff(a, v)
 
-Return the coefficient of `a::HomogeneousPolynomial`, specified by
-`v::Array{Int,1}` which has the indices of the specific monomial.
+Return the coefficient of `a::HomogeneousPolynomial`, specified by `v`,
+which is a tuple (or vector) with the indices of the specific
+monomial.
 """
-function getcoeff(a::HomogeneousPolynomial, v::Array{Int,1})
-    @assert length(v) == get_numvars()
+function getcoeff(a::HomogeneousPolynomial, v::NTuple{N,Int}) where {N}
+    @assert N == get_numvars() && all(v .>= 0)
     kdic = in_base(get_order(),v)
     @inbounds n = pos_table[a.order+1][kdic]
     a[n]
 end
+getcoeff(a::HomogeneousPolynomial, v::Array{Int,1}) = getcoeff(a, (v...,))
 
 getindex(a::HomogeneousPolynomial, n::Int) = a.coeffs[n]
 getindex(a::HomogeneousPolynomial, n::UnitRange) = view(a.coeffs, n)
@@ -123,14 +125,16 @@ setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, c::Colon) where {T<:Number
 """
     getcoeff(a, v)
 
-Return the coefficient of `a::TaylorN`, specified by
-`v::Array{Int,1}` which has the indices of the specific monomial.
+Return the coefficient of `a::TaylorN`, specified by `v`,
+which is a tuple (or vector) with the indices of the specific
+monomial.
 """
-function getcoeff(a::TaylorN, v::Array{Int,1})
+function getcoeff(a::TaylorN, v::NTuple{N, Int}) where {N}
     order = sum(v)
     @assert order ≤ a.order
     getcoeff(a[order], v)
 end
+getcoeff(a::TaylorN, v::Array{Int,1}) = getcoeff(a, (v...,))
 
 getindex(a::TaylorN, n::Int) = a.coeffs[n+1]
 getindex(a::TaylorN, u::UnitRange) = view(a.coeffs, u .+ 1)

src/calculus.jl

@@ -132,10 +132,11 @@ end
 derivative(a::HomogeneousPolynomial, s::Symbol) = derivative(a, lookupvar(s))
 
 """
-    derivative(a, [r=1])
+    derivative(a, r)
 
 Partial differentiation of `a::TaylorN` series with respect
-to the `r`-th variable.
+to the `r`-th variable. The `r`-th variable may be also
+specified through its symbol.
 """
 function derivative(a::TaylorN, r=1::Int)
     T = eltype(a)
@@ -148,6 +149,51 @@ function derivative(a::TaylorN, r=1::Int)
 end
 derivative(a::TaylorN, s::Symbol) = derivative(a, lookupvar(s))
 
+"""
+    derivative(a::TaylorN{T}, ntup::NTuple{N,Int})
+
+Return a `TaylorN` with the partial derivative of `a` defined
+by `ntup::NTuple{N,Int}`, where the first entry is the number
+of derivatives with respect to the first variable, the second is
+the number of derivatives with respect to the second, and so on.
+"""
+function derivative(a::TaylorN, ntup::NTuple{N,Int}) where {N}
+
+    @assert N == get_numvars() && all(ntup .>= 0)
+
+    sum(ntup) > a.order && return zero(a)
+    sum(ntup) == 0 && return copy(a)
+
+    aa = copy(a)
+    for nvar in 1:get_numvars()
+        for numder in 1:ntup[nvar]
+            aa = derivative(aa, nvar)
+        end
+    end
+
+    return aa
+end
+
+"""
+    derivative(ntup::NTuple{N,Int}, a::TaylorN{T})
+
+Returns the value of the coefficient of `a` specified by
+`ntup::NTuple{N,Int}`, multiplied by the corresponding
+factorials.
+"""
+function derivative(ntup::NTuple{N,Int}, a::TaylorN) where {N}
+
+    @assert N == get_numvars() && all(ntup .>= 0)
+
+    c = getcoeff(a, [ntup...])
+    for ind = 1:get_numvars()
+        c *= factorial(ntup[ind])
+    end
+
+    return c
+end
+
+
 ## Gradient, jacobian and hessian
 """
 ```

test/manyvariables.jl

@@ -97,8 +97,8 @@ end
     @test get_order(zeroT) == 1
     @test xT[1][1] == 1
     @test yH[2] == 1
-    @test getcoeff(xT,[1,0]) == 1
-    @test getcoeff(yH,[1,0]) == 0
+    @test getcoeff(xT,(1,0)) == getcoeff(xT,[1,0]) == 1
+    @test getcoeff(yH,(1,0)) == getcoeff(yH,[1,0]) == 0
     @test typeof(convert(HomogeneousPolynomial,1im)) ==
         HomogeneousPolynomial{Complex{Int}}
     @test convert(HomogeneousPolynomial,1im) ==
@@ -154,7 +154,9 @@ end
     @test abs(-1-xT)  == 1+xT
     @test derivative(yH,1) == derivative(xH, :x₂)
     @test derivative(mod2pi(2pi+yT^3),2) == derivative(yT^3, :x₂)
-    @test derivative(yT) == zeroT
+    @test derivative(yT^3, :x₂) == derivative(yT^3, (0,1))
+    @test derivative(yT) == zeroT == derivative(yT, (1,0))
+    @test derivative((0,1), yT) == 1
     @test -xT/3im == im*xT/3
     @test (xH/3im)' == im*xH/3
     @test xT/BigInt(3) == TaylorN(BigFloat,1)/3
@@ -187,11 +189,20 @@ end
     @test integrate(xT^17, :x₁, 2.0) == TaylorN(2.0, get_order())
     @test_throws AssertionError integrate(xT, 1, xT)
     @test_throws AssertionError integrate(xT, :x₁, xT)
-
-
-
-
-    @test derivative(2xT*yT^2,1) == 2yT^2
+    @test_throws AssertionError derivative(xT, (1,))
+    @test_throws AssertionError derivative(xT, (1,2,3))
+    @test_throws AssertionError derivative(xT, (-1,2))
+    @test_throws AssertionError derivative((1,), xT)
+    @test_throws AssertionError derivative((1,2,3), xT)
+    @test_throws AssertionError derivative((-1,2), xT)
+
+
+    @test derivative(2xT*yT^2, (8,8)) == 0
+    @test derivative((8,8), 2xT*yT^2) == 0
+    @test derivative(2xT*yT^2, 1) == 2yT^2
+    @test derivative((1,0), 2xT*yT^2) == 0
+    @test derivative(2xT*yT^2, (1,2)) == 4*one(yT)
+    @test derivative((1,2), 2xT*yT^2) == 4
     @test xT*xT^3 == xT^4
     txy = 1.0 + xT*yT - 0.5*xT^2*yT + (1/3)*xT^3*yT + 0.5*xT^2*yT^2
     @test getindex((1+TaylorN(1))^TaylorN(2),0:4) == txy.coeffs[1:5]
@@ -296,10 +307,10 @@ end
     @test conj(im*yH) == (im*yH)'
     @test conj(im*yT) == (im*yT)'
     @test real( exp(1im * xT)) == cos(xT)
-    @test getcoeff(convert(TaylorN{Rational{Int}},cos(xT)),[4,0]) ==
+    @test getcoeff(convert(TaylorN{Rational{Int}},cos(xT)),(4,0)) ==
         1//factorial(4)
     cr = convert(TaylorN{Rational{Int}},cos(xT))
-    @test getcoeff(cr,[4,0]) == 1//factorial(4)
+    @test getcoeff(cr,(4,0)) == 1//factorial(4)
     @test imag((exp(yT))^(-1im)') == sin(yT)
     exy = exp( xT+yT )
     @test evaluate(exy) == 1
@@ -312,7 +323,7 @@ end
     @test isapprox(evaluate(exy, (1,1)), eeuler^2)
     @test exy(:x₁, 0.0) == exp(yT)
     txy = tan(xT+yT)
-    @test getcoeff(txy,[8,7]) == 929569/99225
+    @test getcoeff(txy,(8,7)) == 929569/99225
     ptxy = xT + yT + (1/3)*( xT^3 + yT^3 ) + xT^2*yT + xT*yT^2
     @test getindex(tan(TaylorN(1)+TaylorN(2)),0:4) == ptxy.coeffs[1:5]
     @test evaluate(xH*yH, 1.0, 2.0) == (xH*yH)(1.0, 2.0) == 2.0
@@ -465,7 +476,12 @@ end
     @test jac == jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] )
     jacobian!(jac, [f1,f2], [2,1])
     @test jac == jacobian([f1,f2], [2,1])
-    @test hessian( f1*f2 ) == [4 -7; -7 0]
+    @test hessian( f1*f2 ) ==
+        [derivative((2,0), f1*f2) derivative((1,1), (f1*f2));
+         derivative((1,1), f1*f2) derivative((0,2), (f1*f2))] == [4 -7; -7 0]
+    @test hessian( f1*f2, [xT, yT] ) ==
+        [derivative(f1*f2, (2,0)) derivative((f1*f2), (1,1));
+         derivative(f1*f2, (1,1)) derivative((f1*f2), (0,2))]
     @test [xT yT]*hessian(f1*f2)*[xT, yT] == [ 2*TaylorN((f1*f2)[2]) ]
     @test hessian(f1^2)/2 == [ [49,0] [0,12] ]
     @test hessian(f1-f2-2*f1*f2) == (hessian(f1-f2-2*f1*f2))'