Remove `Taylor1`? about taylorseries.jl HOT 8 OPEN

I’d like to encourage this proposal.

Making N a type parameter will not make code less efficient: whenever Taylor1 and TaylorN have differing implementations, simply use dispatch so that the most efficient implementation is always used. As structs, their fields are the same, lending themselves readily to parameterisation.

Unifying Taylor1 and TaylorN will simplify a lot of the code in functions.jl especially, large parts of which are inside for T in (:Taylor1, :TaylorN) … end.

And of course, Taylor{N} is more elegant from the user’s perspective, alluding to Array{T,N}.

from taylorseries.jl.

The problem that TaylorN with numvars=1 does not quite mimic Taylor1 behaviour, is related to the way that multivariable Taylor polynomials are constructed and handled. But the example you describe actually works; the subtlety is in the output of set_variables, which is fine in general, but special in the one-variable case, since it returns a vector of one element (and then the error follows). But

julia> x = set_variables("x")  # Returns a vector
1-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0x + 𝒪(‖x‖⁷)

julia> x = set_variables("x")[1]   # let x be the first component
 1.0x₁ + 𝒪(‖x‖⁶)

julia> x*x
 1.0x₁² + 𝒪(‖x‖⁶)

from taylorseries.jl.

Let me now comment on the subtleties of handling a TaylorN object of one variable.

I so think it has a performance cost, because instead of using directly the components of a vector, it has to look up in dictionaries what to use in (trivial) homogeneous polynomials in one variable; yet, I haven't tested this.

But a real subtlety/issue is:

julia> (taylor1_variable(6))^2/taylor1_variable(6)   # this is OK
 1.0⋅t + 𝒪(t⁷)

julia> x^2/x  # x defined as above; should yield x
ERROR: AssertionError: b0 != zero(b0)
 in / at /Users/benet/Fisica/6-IntervalArithmetics/TaylorSeries.jl/src/utils_TaylorN.jl:369

The problem is that I couldn't come with a way to recognize easily that the division of two TaylorN polynomials can be factorized. I can't think (now) of another subtle case like this one.

from taylorseries.jl.

We could probably treat the case Taylor1 as a special case of TaylorN with one variable, as you suggest. But, without thorough tests about the performance hit, I would proceed rather slowly on this.

from taylorseries.jl.

Yes we must obviously check the performance hit.

You are of course right that you have to take the first element of the list if there's only one variable, sorry.
A possibly simpler / better, but also inelegant, solution is

x, = set_variables("x")

A different solution would be to just have a function set_variable for the single variable case.

The division problem occurs also in the following example with 2 variables that one would certainly like to work:

julia> x,y = set_variables("x y")
Warning: redefining constant _params_taylorN
2-element Array{TaylorSeries.TaylorN{Float64},1}:
1.0x + 𝒪(‖x‖⁷)
1.0y + 𝒪(‖x‖⁷)

julia> (x^2_y + x_y^2) / (x*y)
ERROR: AssertionError: b0 != zero(b0)
in / at /Users/dsanders/.julia/v0.4/TaylorSeries/src/utils_TaylorN.jl:369

In this case, the division is by a monomial, so I think it could be checked somehow explicitly if the divisor is a monomial, and then do the division if all of the terms of the Taylor series are divisible by the monomial? This does sound a bit complicated to implement though.

from taylorseries.jl.

I am aware that this kind of things do not work, the reason is the same. I think this is part of the complication that involves representing a polynomial in n variables as vectors of homogeneous polynomials. I was simply unable to solve this. Yet, representing them as a polynomial in one variable, whose coefficients are polynomials in N-1 variables does the trick; this is part of the nice things of #12

EDIT: see this

from taylorseries.jl.

Why not make N a type parameter and then just have a Taylor{N, T} type where N is the number of variables? This way you can specialize on Taylor{1, T} for single variable functions.

I guess this could introduce type instabilities?

from taylorseries.jl.

The short answer is efficiency; Taylor1 is much faster that TaylorNs defined over one variable, as we have them now. It also allows to certain improvements because we are dealing with one variable.

More details now: Though Taylor1 and TaylorN share lots of things, their implementation is a bit different. Taking the product as an example, the Taylor1 case is pretty straight forward; the case for TaylorN is "more or less" the same, if you think that each coefficient of a Taylor1 polynomial is a HomogeneousPolynomial for the multivariable case. Yet, in our current implementation, HomogenousPolynomials are vectors of coefficients, arranged in a specific way that depends on the number of variables. We need to access some tables, for example, to allow to map those locations to the proper monomials, for bookkeeping, and this takes time. There are other limitations related to TaylorN; for example, t^2/t has a well defined Taylor expansion around 0, but x^2/y has no Taylor expansion defined around 0; in other words, it is trickier to factor in the multidimensional case.

Maybe by including the number of the variables as parameters in the struct, as you suggest, Taylor{1,T} can be be made to work as the current Taylor1{T} by dispatch. For me the question is whether we could still have things like Taylor1{TaylorN{T}}, or Taylor{1, Taylor{N,T}}.

from taylorseries.jl.