Matrix evaluation for Taylor1 and TaylorN polynomials

[blas-ko]

Straightforward name! This allows the user to evaluate Taylor-element matrices via evaluate or the function feature for evaluation A(x).

[coveralls]

Coverage decreased (-0.2%) to 97.168% when pulling 6e19f27 on blas-ko:matrix_evaluation into 10ab42c on JuliaDiff:master.

[PerezHz]

Really nice! Works like a charm! 👏

[PerezHz]

travis failure on nightly (macOS) seems to be related to a test which asks for identical equality (test/onevariable.jl line 456):

    @test A*B==Y

The occasional failure of this test is a recurrent issue which has been usually handled by restarting tests, but perhaps it could be fixed if we use an approximate equality, i.e. something like:

    @test norm(A*B-Y,1) ≤ eps(max(norm(A*B,1), norm(Y,1)))
    #or something like: 
    @test norm(A*B-Y,1) ≤ eps()*max(norm(A*B,1), norm(Y,1))

Note that I didn't use isapprox, since A, B and Y are Vector{Taylor1}. So an additional suggestion is to extend isapprox methods for Vector{Taylor1}, etc.

[blas-ko]

Really nice! Works like a charm! 👏

Many thanks! =)

The occasional failure of this test is a recurrent issue which has been usually handled by restarting tests, but perhaps it could be fixed if we use an approximate equality.

Nice that you noticed this! I actually had no idea on why travis was failing!

So an additional suggestion is to extend isapprox methods for Vector{Taylor1}

Maybe this suffices:

import Base: norm, isapprox

norm(v::Vector{T},p::Real=2) where T <: AbstractSeries = norm(norm.(v,p),p)

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

I've tried this with silly examples and does the job, but I'm not sure is the best way to go.

[PerezHz]

@blas-ko do you think a similar thing could be done for 1 and 2 dimensional SubArrays? That would be helpful when working with TaylorIntegration outputs...

[PerezHz]

I've tried this with silly examples and does the job, but I'm not sure is the best way to go.

Well, it LGTM! I think for the moment being you could try that, unless @lbenet or @dpsanders have some other suggestions 😄

[PerezHz]

Upon adding your proposal to the corresponding code files, and substituting @test A*B == Y by either

        @test A*B ≈ Y atol=0.0 rtol=eps()

or

        @test isapprox(A*B, Y, atol=0.0, rtol=eps())

tests pass 👍

[PerezHz]

I mean, tests pass locally

[lbenet]

Thanks for the contribution, LGTM!

I've just restarted the travis job in mac for Julia 0.7 which was failing. Let's see if it passes now and then I'll merge...

[blas-ko]

I've just restarted the travis job in mac for Julia 0.7 which was failing. Let's see if it passes now and then I'll merge...

I'm about to push changes suggested by @PerezHz, so it includes SubArray types in evaluation! Also, should I add isapprox extension in this same PR?

[lbenet]

Ok. Then, I'll wait for the next test round with travis!

[coveralls]

Coverage decreased (-0.1%) to 97.241% when pulling 62bb86b on blas-ko:matrix_evaluation into 10ab42c on JuliaDiff:master.

[lbenet]

Tests are not passing in Julia nightly on Mac, seemingly with an unrelated issue with the matrix multiplication tests; I'm restarting that test...

[lbenet]

By some reason i do not understand, tests fail locally on nighties on a mac with Julia 0.7.0-DEV.3329. The problem appears in matrix multiplication (@test A*B==Y).

Maybe we should allow failures in Julia nightly.

cc: @dpsanders

[blas-ko]

By some reason i do not understand, tests fail locally on nighties on a mac with Julia 0.7.0-DEV.3329

@PerezHz noted that for the A*B == Y test "the occasional failure of this test is a recurrent issue which has been usually handled by restarting tests, but perhaps it could be fixed if we use an approximate equality".

Is it ok if I push the isapprox extension for vectors and see if this solves the inconsistency? @lbenet

[lbenet]

I guess what you mean is to replace line 452 in test/onevariable.jl (which is the problematic test) by something like

@test reduce(&, A*B .≈ Y )

Is that what you mean?

In any case, just go ahead and let's see how travis reacts.

[dpsanders]

Can't you just do evaluate.(A, 3)?

[blas-ko]

Can't you just do evaluate.(A, 3)?

Yes, but it's much slower and it doesn't have the function feature A(3.0)

using TaylorSeries, BenchmarkTools
t = Taylor1(5)
A = [t 2t ; 3t 4t]

@btime evaluate.(A,3.0)
  4.385 μs (24 allocations: 1.19 KiB)

@btime evaluate(A,3.0)
  282.176 ns (12 allocations: 592 bytes)

[dpsanders]

Ah OK, thanks.

Note that you need $ in the call to @btime:

@btime evaluate.($A, 3.0)

[PerezHz]

Thanks @blas-ko for this change! I'd only insist in specifying both atol and rtol as:

@test A*B ≈ Y atol=0.0 rtol=eps()

otherwise, since we are being consistent with Base.isapprox, the default rtol value is sqrt(eps()), which is ~1e-8

[PerezHz]

I guess what you mean is to replace line 452 in test/onevariable.jl (which is the problematic test) by something like @test reduce(&, A*B .≈ Y )

note that @blas-ko defined an isapprox method for vectors of polynomials, so now we can actually do A*B ≈ Y where both the LHS and RHS are vectors of polynomials

[lbenet]

Why not using for the test broadcasting? In current master, the following works:

julia> using TaylorSeries

julia> a = Taylor1(randn(4));

julia> a ≈ a
true

julia> [a] .≈ [a]   # the warnings are related to Julia 0.7
┌ Warning: `rtoldefault(x, y)` is deprecated, use `rtoldefault(x, y, 0)` instead.
│   caller = isapprox at other_functions.jl:189 [inlined]
└                                                  @ Core other_functions.jl:189
1-element BitArray{1}:
 true

julia> [a a] .≈ [a a]
┌ Warning: `rtoldefault(x, y)` is deprecated, use `rtoldefault(x, y, 0)` instead.
│   caller = isapprox at other_functions.jl:189 [inlined]
└                                                  @ Core other_functions.jl:189
1×2 BitArray{2}:
 true  true

The result of .≈ is some sort of array; this is what motivated in my previous comment to suggest @test reduce(&, A*B .≈ Y ), without the need of extending isapprox.

julia> reduce(&, [a a] .≈ [a a])
true

julia> reduce(&, [a a] .≈ [a 2a])
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(uninitialized, m)` instead.
│   caller = *(::Int64, ::Taylor1{Float64}) at arithmetic.jl:215
└                                               @ TaylorSeries arithmetic.jl:215
false

Does extending isapprox make things faster as you showed above for evaluate?

[coveralls]

Coverage decreased (-0.2%) to 97.171% when pulling c22ac64 on blas-ko:matrix_evaluation into 10ab42c on JuliaDiff:master.

[PerezHz]

We could do both:

@test reduce(&, A*B .≈ Y )
@test A*B ≈ Y atol=0.0 rtol=eps()

I agree that for fixing this test we have two solutions and one doesn't involve defining new methods. I just thought that it'd be meaningful to define isapprox for vectors of polynomials

[blas-ko]

I'll push the broadcasting approach on the tests. Do you think is a good idea to extend isapprox as @PerezHz suggested, @lbenet ? I included such extension in my last commit, I can delete that and fix the tests as suggested.

[lbenet]

Let's leave it as it is now, and,if it is not necessary, we'll remove it eventually.

[lbenet]

Merged! Thanks a lot!