JuliaAlgebra / StaticPolynomials.jl

Fast evaluation of multivariate polynomials
https://juliaalgebra.github.io/StaticPolynomials.jl/latest/
Other
16 stars 4 forks source link

Evaluating static polynomials on a Taylor Model #38

Closed ArjunNarayanan closed 4 years ago

ArjunNarayanan commented 4 years ago

I'm trying to use StaticPolynomials.jl along with TaylorModels.jl and ran into an error that I don't quite understand. Example:

import DynamicPolynomials: @polyvar
using StaticPolynomials, TaylorModels

@polyvar x
P = Polynomial(x^2 + 2x + 1)
interval = -1 .. 1
x0 = mid(interval)
order = 1
tm = TaylorModel1(order, x0, interval)

StaticPolynomials.evaluate(P, [tm])

ERROR: LoadError: promotion of types TaylorModel1{Float64,Float64}, Int64 and Int64 failed to change any arguments
Stacktrace:
 [1] sametype_error(::Tuple{TaylorModel1{Float64,Float64},Int64,Int64}) at ./promotion.jl:306
 [2] not_sametype(::Tuple{TaylorModel1{Float64,Float64},Int64,Int64}, ::Tuple{TaylorModel1{Float64,Float64},Int64,Int64}) at ./promotion.jl:300
 [3] promote at ./promotion.jl:289 [inlined]
 [4] muladd(::TaylorModel1{Float64,Float64}, ::Int64, ::Int64) at ./promotion.jl:346
 [5] macro expansion at ./math.jl:100 [inlined]
 [6] macro expansion at /Users/arjun/.julia/packages/StaticPolynomials/1B2wj/src/evaluation.jl:41 [inlined]
 [7] evaluate(::Polynomial{Int64,SExponents{3}(eb76536cb0aa6c4c),Nothing}, ::Array{TaylorModel1{Float64,Float64},1}) at /Users/arjun/.julia/packages/StaticPolynomials/1B2wj/src/evaluation.jl:48
 [8] top-level scope at /Users/arjun/Documents/Berkeley/Research/CodeExamples/FunctionBounds/mwe_static_poly_taylor.jl:11
in expression starting at /Users/arjun/Documents/Berkeley/Research/CodeExamples/FunctionBounds/mwe_static_poly_taylor.jl:11

I'm not sure which is the appropriate place to report this error. Is it an issue with promotion of the TaylorModels type? Any inputs will be helpful!

saschatimme commented 4 years ago

Do you only need to evaluate univariate polynomials? In this case, the @horner macro actually accomplishes the same thing as this package Otherwise this looks like a problem with Taylor Models

ArjunNarayanan commented 4 years ago

I was able to get around this problem by intercepting the muladd(tm::TaylorModelN,a,b) method and overloading it to simply return a*tm+b. I imagine there is a performance hit due to this, but at least it works!

I looked for the documentation of StaticPolynomials.@horner and found StaticPolynomials.@horner_deriv -- is this what you are referring to?

ArjunNarayanan commented 4 years ago

The reason I am using TaylorModels is because I also get a remainder bound for the error. I require this for my application. Does @horner provide this as well?

saschatimme commented 4 years ago

This is in Base, although it's actually called @evalpoly:

help?> Base.@evalpoly
  @evalpoly(z, c...)

  Evaluate the polynomial \sum_k c[k] z^{k-1} for the coefficients c[1], c[2], ...; that is, the coefficients are given in ascending order by power
  of z. This macro expands to efficient inline code that uses either Horner's method or, for complex z, a more efficient Goertzel-like algorithm.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> @evalpoly(3, 1, 0, 1)
  10

  julia> @evalpoly(2, 1, 0, 1)
  5

  julia> @evalpoly(2, 1, 1, 1)
  7

This is the relevan math: https://en.wikipedia.org/wiki/Horner%27s_method

ArjunNarayanan commented 4 years ago

I see, thanks!

To provide some more context -- I am constructing a basis of Lagrange polynomials as a StaticPolynomials.PolynomialSystem. Then for a given set of interpolation coefficients, I am trying to obtain tight bounds for the range of the interpolated polynomial.

Interval arithmetic gives valid bounds, but they are not very tight; hence my interest in Taylor Models.

I guess I can close this issue now since overloading muladd gets around the issue for now!