Add function-like behavior for Taylor1, TaylorN, HomogeneousPolynomial

docs/src/userguide.md

@@ -142,6 +142,23 @@ eBig = evaluate( exp(tBig), one(BigFloat) )
 e - eBig
 ```
 
+Another way to obtain the value of a `Taylor1` polynomial `p` at a given value `x`, is to call `p` as if it were a function:
+
+```@repl userguide
+t = Taylor1(15)
+p = sin(t)
+evaluate(p, pi/2) #get value of p at pi/2 using `evaluate`
+p(pi/2) #get value of p at pi/2 by calling p as a function
+p(pi/2) == evaluate(p, pi/2)
+p(0.0) #get value of `p` at 0.0 by calling p as a function
+evaluate(p) #get p 0-th order coefficient using `evaluate`
+p() #a shortcut to get 0-th order coefficient of `p`
+p() == evaluate(p)
+p() == p(0.0)
+```
+
+Note that the syntax `p(x)` is equivalent to `evaluate(p,x)`, whereas `p()` is equivalent to `evaluate(p)`. For more details about function-like behavior for a given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1) section of the Julia manual.
+
 ## Many variables
 
 A polynomial in ``N>1`` variables can be represented in (at least) two ways:
@@ -283,6 +300,15 @@ number of independent variables.
 evaluate(exy, [.1,.02]) == e^0.12
 ```
 
+Analogously to `Taylor1`, another way to obtain the value of a `TaylorN` polynomial `p` at a given point `x`, is to call `p` as if it were a function:
+
+```@repl userguide
+exy([.1,.02])
+exy([.1,.02]) == e^0.12
+```
+
+Again, note that the syntax `p(x)` for `p::TaylorN` is equivalent to `evaluate(p,x)`, whereas `p()` is equivalent to `evaluate(p)`. For more details about function-like behavior for a given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1) section of the Julia manual.
+
 Functions to compute the gradient, Jacobian and
 Hessian have also been implemented. Using the
 polynomials ``p = x^3 + 2x^2 y - 7 x + 2`` and ``q = y-x^4`` defined above,
@@ -305,7 +331,7 @@ Other specific applications are described in the next [section](examples).
 ## Mixtures
 
 As mentioned above, `Taylor1{T}`, `HomogeneousPolynomial{T}` and `TaylorN{T}`
-are parameterized structures structures such that `T<:AbstractSeries`, the latter
+are parameterized structures such that `T<:AbstractSeries`, the latter
 is a subtype of `Number`. Then, we may actually define Taylor expansions in
 ``N+1`` variables, where one of the variables (the `Taylor1` variable) is
 somewhat special.

src/constructors.jl

@@ -28,6 +28,9 @@ DataType for polynomial expansions in one independent variable.
 - `coeffs :: Array{T,1}` Expansion coefficients; the $i$-th
     component is the coefficient of degree $i-1$ of the expansion.
 - `order  :: Int64` Maximum order (degree) of the polynomial.
+
+Note that `Taylor1` variables are callable. For more information, see
+[`evaluate`](@ref).
 """
 immutable Taylor1{T<:Number} <: AbstractSeries{T}
     coeffs :: Array{T,1}
@@ -81,6 +84,9 @@ DataType for homogenous polynomials in many (>1) independent variables.
 polynomial; the $i$-th component is related to a monomial, where the degrees
 of the independent variables are specified by `coeff_table[order+1][i]`.
 - `order   :: Int` order (degree) of the homogenous polynomial.
+
+Note that `HomogeneousPolynomial` variables are callable. For more information,
+see [`evaluate`](@ref).
 """
 immutable HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}
     coeffs  :: Array{T,1}
@@ -137,6 +143,9 @@ DataType for polynomial expansions in many (>1) independent variables.
 `HomogeneousPolynomial` entries. The $i$-th component corresponds to the
 homogeneous polynomial of degree $i-1$.
 - `order   :: Int`  maximum order of the polynomial expansion.
+
+Note that `TaylorN` variables are callable. For more information, see
+[`evaluate`](@ref).
 """
 immutable TaylorN{T<:Number} <: AbstractSeries{T}
     coeffs  :: Array{HomogeneousPolynomial{T},1}

src/evaluate.jl

@@ -11,7 +11,9 @@
     evaluate(a, [dx])
 
 Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is
-ommitted, its value is considered as zero.
+ommitted, its value is considered as zero. Note that a `Taylor1` polynomial `a`
+may also be evaluated by calling it as a function; that is, the syntax `a(dx)`
+is equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
     @inbounds suma = a[end]
@@ -20,6 +22,7 @@ function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
     end
     suma
 end
+
 function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
     R = promote_type(T,S)
     @inbounds suma = convert(R, a[end])
@@ -28,24 +31,29 @@ function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
     end
     suma
 end
+
 evaluate{T<:Number}(a::Taylor1{T}) = a[1]
 
 doc"""
     evaluate(x, δt)
 
 Evaluates each element of `x::Array{Taylor1{T},1}`, representing
-the dependent variables of an ODE, at *time* δt.
+the dependent variables of an ODE, at *time* δt. Note that an array `x` of
+`Taylor1` polynomials may also be evaluated by calling it as a function;
+that is, the syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()`
+is equivalent to `evaluate(x)`.
 """
 function evaluate{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S)
     R = promote_type(T,S)
     return evaluate(convert(Array{Taylor1{R},1},x), convert(R,δt))
 end
+
 function evaluate{T<:Number}(x::Array{Taylor1{T},1}, δt::T)
     xnew = Array{T}( length(x) )
     evaluate!(x, δt, xnew)
     return xnew
 end
-evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate(a, zero(T))
+evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate.(a)
 
 doc"""
     evaluate!(x, δt, x0)
@@ -62,6 +70,7 @@ function evaluate!{T<:Number}(x::Array{Taylor1{T},1}, δt::T, x0::Union{Array{T,
     end
     nothing
 end
+
 function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Union{Array{T,1},SubArray{T,1}})
     @assert length(x) == length(x0)
     @inbounds for i in eachindex(x)
@@ -70,11 +79,11 @@ function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Uni
     nothing
 end
 
-
 """
     evaluate(a, x)
 
 Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial.
+Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`.
 """
 evaluate{T<:Number,S<:Number}(a::Taylor1{T}, x::Taylor1{S}) =
     evaluate(promote(a,x)...)
@@ -89,7 +98,25 @@ function evaluate{T<:Number}(a::Taylor1{T}, x::Taylor1{T})
     suma
 end
 
+function evaluate{T<:Number}(a::Taylor1{Taylor1{T}}, x::Taylor1{T})
+    @inbounds suma = a[end]
+    @inbounds for k = a.order:-1:1
+        suma = suma*x + a[k]
+    end
+    suma
+end
+
+evaluate{T<:Number,S<:Number}(p::Taylor1{T},x::Array{S}) = evaluate.(p,x)
+
+#function-like behavior for Taylor1
+(p::Taylor1)(x) = evaluate(p, x)
+
+(p::Taylor1)() = evaluate(p)
+
+#function-like behavior for Array{Taylor1,1}
+(p::Array{Taylor1{T},1}){T<:Number}(x) = evaluate(p, x)
 
+(p::Array{Taylor1{T},1}){T<:Number}() = evaluate.(p)
 
 ## Evaluation of multivariable
 function evaluate!{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1},
@@ -100,6 +127,16 @@ function evaluate!{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1},
     end
     nothing
 end
+
+function evaluate!{T<:NumberNotSeriesN}(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1},
+        x0::Array{Taylor1{T},1})
+    @assert length(x) == length(x0)
+    @inbounds for i in eachindex(x)
+        x0[i] = evaluate( x[i], δx )
+    end
+    nothing
+end
+
 function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T,
         x0::Array{TaylorN{T},1})
     @assert length(x) == length(x0)
@@ -110,9 +147,11 @@ function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T,
 end
 
 """
-    evaluate(a, vals)
+    evaluate(a, [vals])
 
-Evaluate a `HomogeneousPolynomial` polynomial at `vals`.
+Evaluate a `HomogeneousPolynomial` polynomial at `vals`. If `vals` is ommitted,
+it's evaluated at zero. Note that the syntax `a(vals)` is equivalent to
+`evaluate(a, vals)`; and `a()` is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T},
         vals::Array{S,1} )
@@ -133,14 +172,21 @@ function evaluate{T<:Number,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T},
 
     return suma
 end
+
 evaluate(a::HomogeneousPolynomial) = zero(a[1])
 
+#function-like behavior for HomogeneousPolynomial
+(p::HomogeneousPolynomial)(x) = evaluate(p, x)
+
+(p::HomogeneousPolynomial)() = evaluate(p)
 
 """
     evaluate(a, [vals])
 
 Evaluate the `TaylorN` polynomial `a` at `vals`.
 If `vals` is ommitted, it's evaluated at zero.
+Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()`
+is equivalent to `evaluate(a)`.
 """
 function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
         vals::Array{S,1})
@@ -165,6 +211,7 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
 
     return sum( sort!(suma, by=abs2) )
 end
+
 function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
         vals::Array{Taylor1{S},1})
     @assert length(vals) == get_numvars()
@@ -190,10 +237,53 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T},
     return suma
 end
 
+function evaluate{T<:NumberNotSeries}(a::TaylorN{Taylor1{T}},
+        vals::Array{Taylor1{T},1})
+    @assert length(vals) == get_numvars()
+
+    num_vars = get_numvars()
+    ct = coeff_table
+    a_length = length(a)
+    ord = maximum( get_order.(vals) )
+    suma = Taylor1(zeros(T, ord))
+
+    for homPol in 1:length(a)
+        sun = zero(Taylor1{T})
+        for (i,a_coeff) in enumerate(a.coeffs[homPol].coeffs)
+            tmp = vals[1]^(ct[homPol][i][1])
+            for n in 2:num_vars
+                tmp *= vals[n]^(ct[homPol][i][n])
+            end
+            suma += a_coeff * tmp
+        end
+    end
+
+    return suma
+end
+
 evaluate{T<:Number}(a::TaylorN{T}) = a[1][1]
 
 function evaluate{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1})
     x0 = Array{T}( length(x) )
     evaluate!( x, δx, x0 )
     return x0
 end
+
+function evaluate{T<:NumberNotSeriesN}(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1})
+    x0 = Array{Taylor1{T}}( length(x) )
+    evaluate!( x, δx, x0 )
+    return x0
+end
+
+evaluate{T<:Number}(x::Array{TaylorN{T},1}) = evaluate.(x)
+
+#function-like behavior for TaylorN
+(p::TaylorN)(x) = evaluate(p, x)
+
+(p::TaylorN)() = evaluate(p)
+
+#function-like behavior for Array{TaylorN,1}
+(p::Array{TaylorN{T},1}){T<:Number}(x) = evaluate(p, x)
+
+(p::Array{TaylorN{T},1}){T<:Number}() = evaluate(p)
+

test/manyvariables.jl

@@ -131,6 +131,12 @@ using Base.Test
     @test (xH/complex(0,BigInt(3)))' ==
         im*HomogeneousPolynomial([BigInt(1),0])/3
     @test evaluate(xH) == zero(eltype(xH))
+    @test xH() == zero(eltype(xH))
+    @test xH([1,1]) == evaluate(xH, [1,1])
+    hp = -5.4xH+6.89yH
+    @test hp([1,1]) == evaluate(hp, [1,1])
+    vr = rand(2)
+    @test hp(vr) == evaluate(hp, vr)
 
     @test derivative(2xT*yT^2,1) == 2yT^2
     @test xT*xT^3 == xT^4
@@ -341,4 +347,34 @@ using Base.Test
     @test_throws AssertionError cos(x)/sin(y)
     @test_throws BoundsError xH[20]
     @test_throws BoundsError xT[20]
+
+    #test function-like behavior for TaylorN
+    @test exy() == 1
+    @test exy([0.1im,0.01im]) == exp(0.11im)
+    @test isapprox(exy([1,1]), e^2)
+    @test sin(asin(xT+yT))([1.0,0.5]) == 1.5
+    @test asin(sin(xT+yT))([1.0,0.5]) == 1.5
+    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(rand(2)) == 1
+    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(zeros(2)) == 1
+    dx = set_variables("x", numvars=4, order=10)
+    P = sin.(dx)
+    v = [1.0,2,3,4]
+    for i in 1:4
+        @test P[i](v) == evaluate(P[i], v)
+    end
+    @test P.(fill(v, 4)) == fill(P(v), 4)
+    F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
+    Q = F(v+dx)
+    @test Q.( fill(v, 4) ) == fill(Q(v), 4)
+    vr = map(x->rand(4), 1:4)
+    @test Q.(vr) == map(x->Q(x), vr)
+    for i in 1:4
+        @test P[i]() == evaluate(P[i])
+        @test Q[i]() == evaluate(Q[i])
+    end
+    @test P() == evaluate.(P)
+    @test P() == evaluate(P)
+    @test Q() == evaluate.(Q)
+    @test Q() == evaluate(Q)
+    @test Q[1:3]() == evaluate(Q[1:3])
 end