Multiplication does not work?

[gaperez64]

TaylorModels.jl/src/arithmetic.jl

Line 291 in e840b99

@assert rnegl_order ≤ get_order()

It seems to me that this assertion is always going to be false? Also, judging by the code just below, the idea is to bound the error due to dropping terms with total degree beyond get_order(). That is useless if we don't allow for the added total order to be larger, right?

TL;DR: Is the assertion needed or in the right direction?

[lbenet]

Thanks for opening ths issue for a clarification.

Yes, it seems that it will always result a false from that assertion. A MWE showing your point follows:

julia> using TaylorModels

julia> x = TaylorModelN(1, 6, IntervalBox(0..0, 0..0), IntervalBox(-1..1, -1..1))
 [1, 1] x₁ + [0, 0]

julia> y = TaylorModelN(2, 6, IntervalBox(0..0, 0..0), IntervalBox(-1..1, -1..1))
 [1, 1] x₂ + [0, 0]

julia> get_order(x)
6

julia> x*y
ERROR: AssertionError: rnegl_order ≤ get_order()
Stacktrace:
 [1] *(a::TaylorModelN{2, Interval{Float64}, Float64}, b::TaylorModelN{2, Interval{Float64}, Float64})
   @ TaylorModels ~/.julia/dev/TaylorModels/src/arithmetic.jl:291
 [2] top-level scope
   @ REPL[6]:1

The trick is, to set-up the internal (Taylor) parameters, so it actually works, either through set_variables or by lowering the order of x and y in the example to 3 (the default parameters are two variables, of maximum order 6, i.e., x = TaylorModelN(1, 3, IntervalBox(0..0, 0..0), IntervalBox(-1..1, -1..1)), and similarly for y):

julia> set_variables("x y", order=12)
2-element Vector{TaylorN{Float64}}:
  1.0 x + 𝒪(‖x‖¹³)
  1.0 y + 𝒪(‖x‖¹³)

julia> x = TaylorModelN(1, 6, IntervalBox(0..0, 0..0), IntervalBox(-1..1, -1..1))
 [1, 1] x + [0, 0]

julia> y = TaylorModelN(2, 6, IntervalBox(0..0, 0..0), IntervalBox(-1..1, -1..1))
 [1, 1] y + [0, 0]

julia> x*y
 [1, 1] x y + [0, 0]

Soooorry, for the lack of documentation. It's lack of time.

[gaperez64]

Thanks!

Yes, I figured changing the internal parameters would make this work. However, it still seems unsatisfactory because the result of multiplication will be truncated at the sum of the total degrees of the operands. If we remove the assertion and have rnegl_order = get_order() as in line 291, then the global order is where we truncate.

Thoughts?

[lbenet]

Mmmmm, not sure I follow you in this...

Taylor models are defined with a polynomial of given order, say n, and the idea is to stick to that order. If two polynomials of that order are multiplied, the result is a polynomial of order 2n. If you impose that your Taylor model resulting from that multiplication has the same order as the factors, you indeed have to truncate to order n, and put the "missing part" into the residual; that's precisely what we do. We follow on that both Berz and Makino as well as Joldes.

If you impose rnegl_order = get_order() you might end computing some trivial zeros, as in the case of x*y I used above. So, one could say that our approach improves performance... But maybe I am too optimistic on this.

[gaperez64]

Yes, I agree this incurs a smaller error.

The alternative (for my purpose) is to have a truncate operation like you do for RTaylorModel1 but I did not find one for TaylorModelN. If it's just because of lack of time, maybe I will try to code one up and send a PR.

[lbenet]

I can't recall now why we have it for the 1-variable cases, and not for the TaylorModelN, but yes, please, if you have time go ahead and submit a PR. That will be certainly very much approceiated!

[lbenet]

Incidentally: a place where we silently change the internal parameters is in the validated integration methods, precisely to account for products and powers (^) which there may appear.

[gaperez64]

Maybe the assertion should just be removed then? I don't see it there for TaylorModel1.