Evaluate for TaylorN fails with ModelingToolkit

[dpsanders]

julia> using TaylorSeries

julia> x, y = set_variables("x y");

julia> t = x^2 + y;

julia> using ModelingToolkit

julia> @variables xx, yy;

julia> t([xx, yy])
ERROR: StackOverflowError:
Stacktrace:
 [1] evaluate(::TaylorN{Float64}, ::Tuple{Operation,Operation}) at /Users/dpsanders/.julia/packages/TaylorSeries/XWMlP/src/evaluate.jl:237 (repeats 79984 times)

[ChrisRackauckas]

Fun

[lbenet]

I do reproduce the problem. No idea where to look for the solution...

[lbenet]

Yet, it does work with HomogeneousPolynomials:

julia> x[1]([xx, yy])
0 + 1.0 * (xx ^ 1 * yy ^ 0)

So the problem may be related to the Horner implementation.

[dpsanders]

I don't think Horner is used for multivariate polynomials?
Shouldn't it just be

sum(h([x, y]) for h in t)` 

to evaluate a TaylorN t?
Or am I missing something?

[dpsanders]

Ah you mean that's too naive and there's some complicated extension of Horner for multivariate polynomials.

[lbenet]

The way we do it is not really Horner, you are right, but some sort of generalization. If I recall correctly, evaluate separately the homogeneous polynomials of given order, sort them, and then add them. Sorting them is to avoid the usual cancellations.

In any case, the observation holds: it seems evaluation works for HomogenousPolynomials.

[lbenet]

I think I found what the problem and it has nothing to do with the way we compute the evaluation

x([xx, yy]) is equivalent to evaluate(x, [xx, yy]) which uses the method defined here, so it calls evaluate(x, (xx, yy)), which uses the same method again and again, until the stack overflow happens:

julia> using ModelingToolkit

julia> using TaylorSeries

julia> x, y = set_variables("x", numvars=2);

julia> @variables xx, yy;

julia> @which evaluate(x, [xx, yy])
evaluate(a::TaylorN, vals) in TaylorSeries at /Users/benet/.julia/dev/TaylorSeries/src/evaluate.jl:237

julia> @which evaluate(x, (xx, yy))
evaluate(a::TaylorN, vals) in TaylorSeries at /Users/benet/.julia/dev/TaylorSeries/src/evaluate.jl:237

Ideally, evaluate(x, (xx, yy)) should use the method at line 222, but it doesn't because S<:NumberNotSeries (which is the type of the NTuple) does not include Expression.

julia> TaylorSeries.NumberNotSeries
Union{Real, Complex}

The question is: how can we update (dynamically) NumberNotSeries so it includes whatever is Number, except for AbstractSeries? NumberNotSeries is defined here.

[dpsanders]

Great catch!

I believe there should not be a NumberNotSeries type at all. Instead there should be

f(::Number)

and

f(::AbstractSeries)

Julia will take care of choosing the correct one by dispatch.

[lbenet]

Using #247 partially solves this issue:

julia> using TaylorSeries, ModelingToolkit

julia> x, y = set_variables("x", numvars=2);

julia> @variables xx, yy;

julia> @which evaluate(x, [xx, yy])
evaluate(a::TaylorN, vals) in TaylorSeries at /Users/benet/.julia/dev/TaylorSeries/src/evaluate.jl:235

julia> @which evaluate(x, (xx, yy))
evaluate(a::TaylorN{T}, vals::Tuple{Vararg{T,N}} where T where N) where T<:Number in TaylorSeries at /Users/benet/.julia/dev/TaylorSeries/src/evaluate.jl:223

Note that the different methods are properly recognized.

Yet, the evaluation doesn't quite work:

julia> evaluate(x, [xx, yy])
ERROR: MethodError: no method matching isless(::Operation, ::Operation)
...

The problem is that evaluating on TaylorNs uses sort! to rearrange the terms before adding them, which I think was introduced for accuracy. We could add a keyword parameter to evaluate for not using it; I prefer this than removing it altogether, because it impacts in the accuracy.

@dpsanders What do you think?

[lbenet]

The following is a proof-of-concept of what I would get using the keyword variable with the example you outlined above (it uses the last commit in #247):

julia> evaluate(t, [xx, yy], sorting=false)
(((((0 + (0 + 1.0 * (xx ^ 0 * yy ^ 1))) + (0 + 1.0 * (xx ^ 2 * yy ^ 0))) + 0) + 0) + 0) + 0

[dpsanders]

That looks great, thanks!

[dpsanders]

What is sorting for?

[lbenet]

I think it was introduced there to add the different contributions of the monomials reordered from the smallest to the largest (according to abs2)... The old story of including properly as many as possible significant figures.

[lbenet]

Reopening, as a reminder...