Missing functionality: converting coefficients

[nsajko]

I have some code where I define multiple AbstractPolynomial subtypes. In some cases I'd like to be able to turn a polynomial over one ring to a polynomial over another ring, by converting all its coefficients. It would seem appropriate for a function with that functionality to live in the Polynomials module. Something like this:

over_ring(::Type{New}, p::AbstractPolynomial{Old, x})::AbstractPolynomial{New, x} where {New <: Number, Old <: Number, x}

For example, for Polynomial a possible implementation would be this:

over_ring(::Type{New}, p::Polynomial{Old, x}) where {New <: Number, Old <: Number, x} =
  Polynomial{New, x}(p)

Usage example:

julia> over_ring(Float64, Polynomial([3, 2, 1]))
Polynomial(3.0 + 2.0*x + 1.0*x^2)

I'd still need to define a method of this function for each subtype of AbstractPolynomial owned by me, but it would be much nicer and tidier if I could overload a Polynomials function.

Would an interface like this would fit well into Polynomials?

[nsajko]

If this is added, it might make sense to also add something analogous for changing variables:

with_variable(::Val{y}, p::Polynomial{R, x})::Polynomial{R, y} where {x, y, R <: Number}

[nsajko]

BTW, at first I thought that Base.similar could be used for this purpose, but then I realized that Polynomials already overloads Base.similar for some other purpose.

Happy New Year!

[jverzani]

Hi

Usage example:

julia> over_ring(Float64, Polynomial([3, 2, 1]))
Polynomial(3.0 + 2.0*x + 1.0*x^2)

I think this is convert:

julia> convert(Polynomial{Float64}, Polynomial([3, 2, 1]))
Polynomial(3.0 + 2.0*x + 1.0*x^2)

But convert fails when a different symbol is attempted. I don't think changing that is a great idea, as I'm not sure it would be widely used. A workaround might compose coeffs with the constructor:

julia> Polynomial{Float64, :u}(coeffs(Polynomial([3, 2, 1])))
Polynomial(3.0 + 2.0*u + 1.0*u^2)

or convert could also be used:

julia> convert(Polynomial{Float64,:u}, coeffs(Polynomial([3, 2, 1])))
Polynomial(3.0 + 2.0*u + 1.0*u^2)

Let me know if I didn't understand the use case.

[nsajko]

The crucial difference is that convert takes the polynomial type as an argument, while over_ring would take just the ring. To use convert I need to know the full concrete type of the AbstractPolynomial (the return type), which inhibits generic code.

In this way over_ring for AbstractPolynomial types would be analogous to Base.similar for AbstractArray types.

Suppose you have some function like this: func(p::AbstractPolynomial{Coef}) and suppose we want to convert the coefficient type of p somewhere in the function body. It's not possible to use convert on p because we don't know the concrete type of p stripped of its type parameters, and it could be any subtype of AbstractPolynomial.

In my code the example is like so: I have some type CompressedPolynomial <: AbstractPolynomial:

https://gitlab.com/nsajko/PolynomialApproximation.jl/-/blob/main/src/CompressedPolynomials.jl#L12-L36

I want to be able to change the coefficient type of an instance of CompressedPolynomial. The problem is that each CompressedPolynomial is a wrapper over another type P <: AbstractPolynomial:

https://gitlab.com/nsajko/PolynomialApproximation.jl/-/blob/main/src/CompressedPolynomials.jl#L27

So to solve this using convert, I'd have to special-case each possible wrapped P <: AbstractPolynomial: this would necessitate an explosion of lines of code (a method for each possible AbstractPolynomial). It would also mean I'd have to whitelist supported types for P, instead of accepting any AbstractPolynomial.

So, to solve this generically, I introduced a new module called OverRing:

https://gitlab.com/nsajko/PolynomialApproximation.jl/-/blob/main/src/OverRing.jl#L5-L14

It hosts just the function over_ring, with methods for Polynomial and ImmutablePolynomial. But now I can overload OverRing.over_ring for any of my subtypes of AbstractPolynomial, which enables this trivial implementation of what I wanted in the first place:

https://gitlab.com/nsajko/PolynomialApproximation.jl/-/blob/main/src/CompressedPolynomials.jl#L53

The important part of the line is this expression: over_ring(New, p.main_factor). Even though I know within that method body that the concrete type of p.main_factor is exactly P, there's no way within Julia to strip the type parameters from P, preventing me from using convert instead of over_ring:

https://docs.julialang.org/en/v1/manual/methods/#Building-a-similar-type-with-a-different-type-parameter

There is no general transform of one subtype into another subtype with a different parameter.

So in the end I figured out a way to solve this, because now my code supports any AbstractPolynomial type that additionally implements OverRing.over_ring. But still, a better place for the over_ring function to live would be within Polynomials. This would prevent other people from solving the same problem I had to deal with, and it would also be good for compatibility within the ecosystem if there was a canonical function to overload in cases like these.

I hope I was more clear this time?

[jverzani]

Sorry, you were clear. I just didn't read carefully. I'm thinking of Julia generics for this, and I don't think similar is correct, but reinterpret seems to be, though I doubt you would want a view of the original data, so not quite right. I'll put in a PR with the fix, as it does seem like a good inclusion and we can sort out the naming.

[nsajko]

Two issues with the current Polynomials.copy_with_eltype implementation:

1. The type inference still doesn't work perfectly for ImmutablePolynomial. There's also lots of run time dispatch in that case.
2. I'm pretty sure that the first argument (the coefficient type), unlike the second argument (the variable name), doesn't need to be wrapped in Val. The point of using Val is usually to move values into the type domain, but types are obviously already in the type domain, so using Val instead of Type just makes the usage more cumbersome.

Example for the type inference and run time dispatch issue:

  | | |_| | | | (_| |  |  Version 1.10.0-DEV.982 (2023-04-10)
 _/ |\__'_|_|_|\__'_|  |  Commit ea72b942792 (4 days old master)
|__/                   |

julia> using JET, Test, Polynomials

julia> p = ImmutablePolynomial((3, 2, 1))
ImmutablePolynomial(3 + 2*x + x^2)

julia> Polynomials.copy_with_eltype(Val(Float64), p)
ImmutablePolynomial(3.0 + 2.0*x + 1.0*x^2)

julia> @inferred Polynomials.copy_with_eltype(Val(Float64), p)
ERROR: return type ImmutablePolynomial{Float64, :x, 3} does not match inferred return type ImmutablePolynomial{Float64, :x}
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] top-level scope
   @ REPL[4]:1

julia> @report_opt Polynomials.copy_with_eltype(Val(Float64), p)
═════ 48 possible errors found ═════
┌ @ /home/nsajko/.julia/packages/Polynomials/Fh8md/src/common.jl:539 Polynomials.copy_with_eltype(V, Polynomials.Val(Y), p)
[...]

[jverzani]

Thanks. For 2., I don't have a good argument as to why I wrapped the T in Val. For 1. I should double check, by my guess is that N parameter is the problem as not all operations keep that the same. Let me play around with Jet, a tool I haven't taken advantage of yet, and see if there is low hanging fruit to be found.

[nsajko]

Yeah, your guess is right. To fix (1), I think you could simply special-case the implementation for ImmutablePolynomial, similarly as shown in my link above from January:

https://gitlab.com/nsajko/PolynomialApproximation.jl/-/blob/main/src/OverRing.jl#L12-L14

[jverzani]

I think #492 is a better take on this issue. Thanks for your help.