Better compatibility with packages such as ApproxFun.jl and Polynomials.jl

[benjione]

If DynamicPolynomials library is used together with e.g. ApproxFun, there is an error if the package IntervalSets is used and the following is executed with v being of type PolyVar:

_in(v, I::TypedEndpointsInterval{:closed,:closed}) = leftendpoint(I) ≤ v ≤ rightendpoint(I)
_in(v, I::TypedEndpointsInterval{:open,:open}) = leftendpoint(I) < v < rightendpoint(I)
_in(v, I::TypedEndpointsInterval{:closed,:open}) = leftendpoint(I) ≤ v < rightendpoint(I)
_in(v, I::TypedEndpointsInterval{:open,:closed}) = leftendpoint(I) < v ≤ rightendpoint(I)

The error is because Base.isless is not defined for PolyVar
Code example producing the error:

using SumOfSquares
using DynamicPolynomials
using DynamicPolynomials: PolyVar
using ApproxFun
using SCS
using IntervalSets

## use it for something:
@polyvar x

model = SOSModel(SCS.Optimizer)
@variable(model, a[1:10])
p = Fun(Chebyshev(), a)
p(x)  # this throws for example the error

One can manually overwrite those functions for the desired type, e.g.

using DynamicPolynomials: PolyVar
using IntervalSets

## define these functions as a trick:
IntervalSets.in(v::PolyVar, I::IntervalSets.TypedEndpointsInterval{:closed,:closed}) = true
IntervalSets.in(v::PolyVar, I::IntervalSets.TypedEndpointsInterval{:open,:open}) = true
IntervalSets.in(v::PolyVar, I::IntervalSets.TypedEndpointsInterval{:closed,:open}) = true
IntervalSets.in(v::PolyVar, I::IntervalSets.TypedEndpointsInterval{:open,:closed}) = true

Is it possible to overwrite the IntervalSets.in in order to get better compatibility between the polynomial packages?
Or has this to be done in MultivariatePolynomials.jl instead of here?

[benjione]

Also see here and this pull request here for how to do this specialization.

[blegat]

I'm tempted to say that we don't want to define in for the same reason we do not define isless. It does not make sense to ask whether a symbolic variable is less than a number. Maybe there is a function higher in the stacktrace that would make sense to define

[benjione]

While this is true, I think the desired behavior would be:

julia> @polyvar x
(x,)
julia> 0 ≤ x ≤ 1 # this would have a side effect on x
0 ≤ x ≤ 1
julia> p = Fun(Chebyshev(), [0.0, 1.0, 2.0])
Fun(Chebyshev(), [0.0, 1.0, 2.0])
juliy> p(x)
4x^2 + x - 2 for 0 ≤ x ≤ 1

Or maybe something like this to avoid side effects:

julia> @polyvar x
(x,)
julia> x2 = impose(x, lowerbound=0, upperbound=1) # no side effect on x, Base.isless can be defined as x2 < y if true for all x2
0 ≤ x2 ≤ 1
julia> @impose 0 ≤ x ≤ 1 # maybe a macro overwriting x
0 ≤ x ≤ 1
julia> p = Fun(Chebyshev(), [0.0, 1.0, 2.0])
Fun(Chebyshev(), [0.0, 1.0, 2.0])
juliy> p(x2)
4x2^2 + x2 - 2 for 0 ≤ x2 ≤ 1

[blegat]

We have such concept of domain for x in https://github.com/JuliaAlgebra/SemialgebraicSets.jl but it wouldn't work like that.
Note that this can be achieved as follows:

julia> @polyvar x
(x,)
julia> p = Fun(Chebyshev(), [0.0, 1.0, 2.0])
Fun(Chebyshev(), [0.0, 1.0, 2.0])
julia> clenshaw(p.space, p.coefficients, x)
4x2^2 + x2 - 2

See https://github.com/JuliaApproximation/ApproxFunOrthogonalPolynomials.jl/blob/8bafeff745ad5191916cf26d909a152cd29b7c48/src/Spaces/PolynomialSpace.jl#L21
That's a bit hacky though. @jishnub is there a recommended way in ApproxFun to evaluate while skipping the domain check ?

[jishnub]

I don't think there's a way at present, although this does sound like a good idea. Perhaps @benjione can do this at present (an extension of their impose idea):

julia> struct InDomain{T}
       x :: T
       end

julia> p = Fun(Chebyshev(), [0.0, 1.0, 2.0])
Fun(Chebyshev(), [0.0, 1.0, 2.0])

julia> ApproxFunBase.tocanonical(::ChebyshevInterval{<:Number}, x::InDomain) = x.x

julia> Base.in(::InDomain, ::ChebyshevInterval) = true

julia> p(InDomain(x))
4.0x² + x - 2.0

However, using clenshaw directly should be fine as well:

julia> clenshaw(space(p), ApproxFun.coefficients(p), x)
4.0x² + x - 2.0

[benjione]

Thanks, I think I will use the InDomain{T} solution for now.

Nevertheless, I think it would be nice if polyvar would work with other polynomial packages. E.g. Polynomials.jl throws exact the same error as ApproxFun.jl.

Maybe the InDomain{T} solution could be added somewhere in the documentation.
But probably there are not many users anyway combining these packages.

[jishnub]

E.g. Polynomials.jl throws exact the same error as ApproxFun.jl.

what's the error? This seems to work for me:

julia> @polyvar x
(x,)

julia> Polynomials.Polynomial([1,2,3])(x)
3x² + 2x + 1

[benjione]

I think the normal polynomials are not domain bounded. For Chebyshev, I get:

ChebyshevT(rand(10))(x)
ERROR: MethodError: no method matching isless(::Int64, ::PolyVar{true})
Closest candidates are:
...

[blegat]

We try to avoid defining methods just to make a workflow work as it might break another workflow. Here, defining isless between a variable and a polynomial variable would be confusing for other workflow so we'd rather not do it.

[dlfivefifty]

I'm adding support for ClassicalOrthogonalPolynomials.jl (which is meant to underpin ApproxFun.jl in the future.... scheduled for 2050s 😅):

JuliaApproximation/ClassicalOrthogonalPolynomials.jl#214

[dlfivefifty]

Btw, is there any support for Rational coefficients? I don't think I need it but usually algebraicists prefer it ...

[blegat]

Thanks for the new package, this could be useful for https://github.com/JuliaAlgebra/MultivariateBases.jl
Yes, the coefficients can be anything, it does not even have to be a subtype of Number, I make sure it even works when the coefficients are JuMP expressions or MOI functions as this is needed by SumOfSquares.jl

[dlfivefifty]

OK I was expecting this to work but maybe I should file an issue:

julia> using DynamicPolynomials

julia> @polyvar x
(x,)

julia> x/2
0.5x

julia> x//2
ERROR: MethodError: no method matching //(::Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, ::Int64)
The function `//` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  //(::Integer, ::Integer)
   @ Base rational.jl:84
  //(::Rational, ::Integer)
   @ Base rational.jl:86
  //(::Complex, ::Real)
   @ Base rational.jl:100
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[4]:1