Matrix evaluation for Taylor1 and TaylorN polynomials (#144)

* Added matrix evaluation for Taylor1 and TaylorN.

* Added tests for matrix evaluation.

* Extended evaluation methods to include SubArrays.

* Deleted redundant functions.

* Updated Matrix evaluation tests.

* extended norm and isapprox for AbstractSeries vectors.

* Fixed problematic test

src/evaluate.jl

@@ -35,13 +35,13 @@ evaluate(a::Taylor1{T}) where {T<:Number} = a[0]
 doc"""
     evaluate(x, δt)
 
-Evaluates each element of `x::Array{Taylor1{T},1}`, representing
+Evaluates each element of `x::Union{ Vector{Taylor1{T}}, Matrix{Taylor1{T}} }`, representing
 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(x::Array{Taylor1{T},1}, δt::S) where {T<:Number, S<:Number}
+function evaluate(x::Union{Array{Taylor1{T},1}, SubArray{Taylor1{T},1}}, δt::S) where {T<:Number, S<:Number}
     R = promote_type(T,S)
     return evaluate(convert(Array{Taylor1{R},1},x), convert(R,δt))
 end
@@ -50,7 +50,28 @@ function evaluate(x::Array{Taylor1{T},1}, δt::T) where {T<:Number}
     evaluate!(x, δt, xnew)
     return xnew
 end
+
 evaluate(a::Array{Taylor1{T},1}) where {T<:Number} = evaluate(a, zero(T))
+evaluate(a::SubArray{Taylor1{T},1}) where {T<:Number} = evaluate(a, zero(T))
+
+function evaluate(A::Union{Array{Taylor1{T},2}, SubArray{Taylor1{T},2}}, δt::S) where {T<:Number, S<:Number}
+    R = promote_type(T,S)
+    return evaluate(convert(Array{Taylor1{R},2},A), convert(R,δt))
+end
+function evaluate(A::Array{Taylor1{T},2}, δt::T) where {T<:Number}
+    n,m = size(A)
+    Anew = Array{T}( n,m )
+    xnew = Array{T}( n )
+
+    for i in 1:m
+        evaluate!(A[:,i], δt, xnew)
+        Anew[:,i] = xnew
+    end
+
+    return Anew
+end
+evaluate(A::Array{Taylor1{T},2}) where {T<:Number} = evaluate.(A)
+evaluate(A::SubArray{Taylor1{T},2}) where {T<:Number} = evaluate.(A)
 
 doc"""
     evaluate!(x, δt, x0)
@@ -115,10 +136,17 @@ evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} =
 
 (p::Taylor1)() = evaluate(p)
 
-#function-like behavior for Array{Taylor1,1}
+#function-like behavior for Vector{Taylor1}
 (p::Array{Taylor1{T},1})(x) where {T<:Number} = evaluate(p, x)
-
+(p::SubArray{Taylor1{T},1})(x) where {T<:Number} = evaluate(p, x)
 (p::Array{Taylor1{T},1})() where {T<:Number} = evaluate.(p)
+(p::SubArray{Taylor1{T},1})() where {T<:Number} = evaluate.(p)
+
+#function-like behavior for Matrix{Taylor1}
+(p::Array{Taylor1{T},2})(x) where {T<:Number} = evaluate(p, x)
+(p::SubArray{Taylor1{T},2})(x) where {T<:Number} = evaluate(p, x)
+(p::Array{Taylor1{T},2})() where {T<:Number} = evaluate.(p)
+(p::SubArray{Taylor1{T},2})() where {T<:Number} = evaluate.(p)
 
 ## Evaluation of multivariable
 function evaluate!(x::Array{TaylorN{T},1}, δx::Array{T,1},
@@ -304,36 +332,53 @@ end
 
 evaluate(a::TaylorN{T}) where {T<:Number} = a[0][1]
 
+#Vector evaluation
+function evaluate(x::Union{Array{TaylorN{T},1},SubArray{TaylorN{T},1}}, δx::Vector{S}) where {T<:Number, S<:Number}
+    R = promote_type(T,S)
+    return evaluate(convert(Array{TaylorN{R},1},x), convert(Vector{R},δx))
+end
+
 function evaluate(x::Array{TaylorN{T},1}, δx::Array{T,1}) where {T<:Number}
     x0 = Array{T}( length(x) )
     evaluate!( x, δx, x0 )
     return x0
 end
 
-function evaluate(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1}) where
-        {T<:NumberNotSeriesN}
+evaluate(x::Array{TaylorN{T},1}) where {T<:Number} = evaluate.(x)
+evaluate(x::SubArray{TaylorN{T},1}) where {T<:Number} = evaluate.(x)
 
-    x0 = Array{Taylor1{T}}( length(x) )
-    evaluate!( x, δx, x0 )
-    return x0
+#Matrix evaluation
+function evaluate(A::Union{Array{TaylorN{T},2}, SubArray{TaylorN{T},2}}, δx::Vector{S}) where {T<:Number, S<:Number}
+    R = promote_type(T,S)
+    return evaluate(convert(Array{TaylorN{R},2},A), convert(Vector{R},δx))
 end
+function evaluate(A::Array{TaylorN{T},2}, δx::Vector{T}) where {T<:Number}
+    n,m = size(A)
+    Anew = Array{T}( n,m )
+    xnew = Array{T}( n )
+
+    for i in 1:m
+        evaluate!(A[:,i], δx, xnew)
+        Anew[:,i] = xnew
+    end
 
-function evaluate(x::Array{TaylorN{T},1}, δx::Array{TaylorN{T},1}) where
-        {T<:NumberNotSeriesN}
-
-    x0 = Array{TaylorN{T}}( length(x) )
-    evaluate!( x, δx, x0 )
-    return x0
+    return Anew
 end
-
-evaluate(x::Array{TaylorN{T},1}) where {T<:Number} = evaluate.(x)
+evaluate(A::Array{TaylorN{T},2}) where {T<:Number} = evaluate.(A)
+evaluate(A::SubArray{TaylorN{T},2}) where {T<:Number} = evaluate.(A)
 
 #function-like behavior for TaylorN
 (p::TaylorN)(x) = evaluate(p, x)
-
 (p::TaylorN)() = evaluate(p)
 
-#function-like behavior for Array{TaylorN,1}
+#function-like behavior for Vector{TaylorN}
 (p::Array{TaylorN{T},1})(x) where {T<:Number} = evaluate(p, x)
-
+(p::SubArray{TaylorN{T},1})(x) where {T<:Number} = evaluate(p, x)
 (p::Array{TaylorN{T},1})() where {T<:Number} = evaluate(p)
+(p::SubArray{TaylorN{T},1})() where {T<:Number} = evaluate(p)
+
+#function-like behavior for Matrix{TaylorN}
+(p::Array{TaylorN{T},2})(x) where {T<:Number} = evaluate(p, x)
+(p::SubArray{TaylorN{T},2})(x) where {T<:Number} = evaluate(p, x)
+(p::Array{TaylorN{T},2})() where {T<:Number} = evaluate.(p)
+(p::SubArray{TaylorN{T},2})() where {T<:Number} = evaluate.(p)

src/other_functions.jl

@@ -159,6 +159,8 @@ which returns a non-negative number.
 norm(x::AbstractSeries, p::Real=2) = norm( norm.(x.coeffs, p), p)
 norm(x::Union{Taylor1{T},HomogeneousPolynomial{T}}, p::Real=2) where
     {T<:NumberNotSeries} = norm(x.coeffs, p)
+#norm for Taylor vectors
+norm(v::Vector{T}, p::Real=2) where T<:AbstractSeries = norm( norm.(v, p), p)
 
 # rtoldefault
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
@@ -190,7 +192,13 @@ function isapprox(x::T, y::S; rtol::Real=rtoldefault(x,y), atol::Real=0.0,
         norm(x-y,1) <= atol + rtol*max(norm(x,1), norm(y,1))) ||
         (nans && isnan(x) && isnan(y))
 end
+#isapprox for vectors of Taylors
+function isapprox(x::Vector{T}, y::Vector{S}; rtol::Real=rtoldefault(T,S), atol::Real=0.0,
+        nans::Bool=false) where {T<:AbstractSeries, S<:AbstractSeries}
 
+    x == y || norm(x-y,1) <= atol + rtol*max(norm(x,1), norm(y,1)) ||
+        (nans && isnan(x) && isnan(y))
+end
 
 #taylor_expand function for Taylor1
 doc"""
@@ -267,7 +275,7 @@ end
 for T in (:Taylor1, :TaylorN)
     @eval deg2rad(z::$T{T}) where {T<:AbstractFloat} = z * (convert(T, pi) / 180)
     @eval deg2rad(z::$T{T}) where {T<:Real} = z * (convert(float(T), pi) / 180)
-    
+
     @eval rad2deg(z::$T{T}) where {T<:AbstractFloat} = z * (180 / convert(T, pi))
     @eval rad2deg(z::$T{T}) where {T<:Real} = z * (180 / convert(float(T), pi))
 end

test/manyvariables.jl

@@ -279,6 +279,13 @@ end
     v = zeros(Int, 2)
     @test evaluate!([xT, yT], ones(Int, 2), v) == nothing
     @test v == ones(2)
+    A_TN = [xT 2xT 3xT; yT 2yT 3yT]
+    @test evaluate(A_TN, ones(2)) == [1.0 2.0 3.0; 1.0 2.0 3.0]
+    @test A_TN() == [0.0  0.0  0.0; 0.0  0.0  0.0]
+    t = Taylor1(10)
+    @test A_TN([t,t^2]) == [t 2t 3t; t^2 2t^2 3t^2]
+    @test view(A_TN, :, :)(ones(2)) == A_TN(ones(2))
+    @test view(A_TN, :, 1)(ones(2)) == A_TN[:,1](ones(2))
 
     @test evaluate(sin(asin(xT+yT)), [1.0,0.5]) == 1.5
     @test evaluate(asin(sin(xT+yT)), [1.0,0.5]) == 1.5

test/onevariable.jl

@@ -230,6 +230,13 @@ end
     taylor_x = exp(Taylor1(Float64,13))
     @which evaluate(taylor_x, taylor_a)
     @test taylor_x(taylor_a) == evaluate(taylor_x, taylor_a)
+    A_T1 = [t 2t 3t; 4t 5t 6t ]
+    @test evaluate(A_T1,1.0) == [1.0  2.0  3.0; 4.0  5.0  6.0]
+    @test evaluate(A_T1,1.0) == A_T1(1.0)
+    @test evaluate(A_T1) == A_T1()
+    @test A_T1(tsquare) == [tsquare 2tsquare 3tsquare; 4tsquare 5tsquare 6tsquare]
+    @test view(A_T1, :, :)(1.0) == A_T1(1.0)
+    @test view(A_T1, :, 1)(1.0) == A_T1[:,1](1.0)
 
     @test sin(asin(tsquare)) == tsquare
     @test tan(atan(tsquare)) == tsquare
@@ -449,7 +456,7 @@ end
         A_mul_B!(Y,A,B)
 
         # do we get the same result when using the `A*B` form?
-        @test A*B==Y
+        @test A*B≈Y
         # Y should be extended after the multilpication
         @test reduce(&, [y1.order for y1 in Y] .== Y[1].order)
         # B should be unchanged