Mixing of variables in concatenation operations allowed

andreasvarga commented 3 years ago

In the current version v2.0 the mixing of variables is allowed in the concatenation operations:

julia> x = Polynomial([1, 1])
Polynomial(1 + x)

julia> y = Polynomial([1, 1],:y)
Polynomial(1 + y)

julia> [x y]
1×2 Array{Polynomial{Int64,X} where X,2}:
 Polynomial(1 + x)  Polynomial(1 + y)

Is this a feature or a bug?

I hope very much this is a bug!

jverzani commented 3 years ago

Bug. Sorry, let's see where that needs to be nipped.

jverzani commented 3 years ago

I might have a fix, but let me ask what you would like these to return

[p q] # Error, right?

[p one(q)] # same as [p 1] with symbol of p, right?

[one(p) q] # same as [1 q] with symbol of q, right?

But if we do the last two, then what is?

[one(p) one(q)]

My "fix" would make order dependent as the symbol would come from the first constant polynomial encountered. This also occurs with one(p) + one(q) or one(p) * one(q). Does that seem odd to you?

andreasvarga commented 3 years ago

It is OK to be done in this way. Thus constants are tolerated even they have another symbol.

john verzani @.***> schrieb am Fr., 12. März 2021, 18:29:

I might have a fix, but let me ask what you would like these to return

[p q] # Error, right?

[p one(q)] # same as [p 1] with symbol of p, right?

[one(p) q] # same as [1 q] with symbol of q, right?

But if we do the last two, then what is?

[one(p) one(q)]

My "fix" would make order dependent as the symbol would come from the first constant polynomial encountered. This also occurs with one(p) + one(q) or one(p) * one(q). Does that seem odd to you?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/312#issuecomment-797641027, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEF3OFE4IE3HNDZH2O3TDI6OXANCNFSM4ZCNVYTQ .

jverzani commented 3 years ago

The basic fix is a generalization of this:

Base.promote_eltypeof(p::P, q::Q) where {T,X, P<: Polynomial{T,X}, S, Q<:Polynomial{S,X}} =
    Polynomial{promote_type(T,S),X}

function Base.promote_eltypeof(p::P, q::Q) where {T,X, P<: Polynomial{T,X}, S,Y,Q<:Polynomial{S,Y}} 
    Polynomials.isconstant(p) && return Polynomial{promote_type(T,S), Y}
    Polynomials.isconstant(q) && return Polynomial{promote_type(T,S), X}
    error("X Y don't match")
end

Rather than ask you to go through github to try a branch, could you just execute the above and see if things work as expected for you?

jverzani commented 3 years ago

Not sure this is a rabbit hole to keep chasing, but I find another method seems needed based on this experiment you should be able to try at home:

using Polynomials
using LinearAlgebra
using Test

p = Polynomial([1, 2], :x)
q = Polynomial([1, 2], :y)

[p q] # should error, instead 1×2 Array{Polynomial{Int64,X} where X,2}:
[p one(q)] # should be Array{Ppolynomial{Int64, :x}} is 1×2 Array{Polynomial{Int64,X} where X,2}:

Base.promote_eltypeof(p::P, q::Q) where {T,X, P<: Polynomial{T,X}, S, Q<:Polynomial{S,X}} =
    ⟒(P){promote_type(T,S),X}
function Base.promote_eltypeof(p::P, q::Q) where {T,X, P<:Polynomial{T,X}, S,Y,Q<:Polynomial{S,Y}} 
    Polynomials.isconstant(p) && return Polynomial{promote_type(T,S), Y}
    Polynomials.isconstant(q) && return Polynomial{promote_type(T,S), X}
    Polynomials.assert_same_variable(X,Y)
end

@test_throws ArgumentError [p q]
@test eltype([p one(q)]) == typeof(p)

# various constructs
[p p] # okay Array{Polynomial{Int64,:x},2}
[p; p] # okay

[p 1] # okay
[p one(q)] # okay
[one(p) one(q)] # kinda okay, has :y as a symbol

A = [p p^2; p^3 one(q)] # okay

B = [[p p^2] I] # okay

C = [A I] # okay

D = [p I] # not okay as I does not promote to Polynomial

E = [[p p]; [1 1]] # okay

F = [[p p]; [one(p) one(p)]] # okay

G = [[p p]; [one(p) one(q)]] # not okay; [one(p) one(q)] has symbol :y so type if off

H = [[p p]; [q q]] # not okay Array{Array{T,2} where T,1}:

function Base.promote_eltype(p::PP, q::QQ) where {T,X,N, P<:Polynomial{T,X}, PP<:Array{P,N},
                                                  S,Y,   Q<:Polynomial{S,Y}, QQ<:Array{Q,N}} 
    Polynomials.assert_same_variable(X,Y)
    ⟒(P){promote_type(T,S),X}
end

H = [[p p]; [q q]] # not okay Array{Array{T,2} where T,1}:

@test_throws ArgumentError [[p p]; [q q]]

@test_throws ArgumentError  [[p p]; [one(q) one(q)]]

[[p p] [ 1 one(p)]]

Are there other combinations I should also be looking at?

andreasvarga commented 3 years ago

I found the following errors:

julia> [p q 1]
1×3 Array{Polynomial{Int64,X} where X,2}:
 Polynomial(1 + 2*x)  Polynomial(1 + 2*y)  Polynomial(1)

julia>  [p; q; 1]
3-element Array{Polynomial{Int64,X} where X,1}:
 Polynomial(1 + 2*x)
 Polynomial(1 + 2*y)
 Polynomial(1)

julia> G = [[p p]; [one(p) one(q)]]
2×2 Array{Polynomial{Int64,X} where X,2}:
 Polynomial(1 + 2*x)  Polynomial(1 + 2*x)
 Polynomial(1)        Polynomial(1)

julia> H = [[p p]; [q q]]
2×2 Array{Polynomial{Int64,X} where X,2}:
 Polynomial(1 + 2*x)  Polynomial(1 + 2*x)
 Polynomial(1 + 2*y)  Polynomial(1 + 2*y)

After adding

function Base.promote_eltype(p::PP, q::QQ) where {T,X,N, P<:Polynomial{T,X}, PP<:Array{P,N},
                                                  S,Y,   Q<:Polynomial{S,Y}, QQ<:Array{Q,N}} 
    Polynomials.assert_same_variable(X,Y)
    ⟒(P){promote_type(T,S),X}
end

I got


julia> G = [[p p]; [one(p) one(q)]]
ERROR: ArgumentError: Polynomials have different indeterminates
Stacktrace:
 [1] assert_same_variable(::Symbol, ::Symbol) at C:\Users\Andreas\.julia\packages\Polynomials\1nhQ0\src\common.jl:394
 [2] promote_eltype(::Array{Polynomial{Int64,:x},2}, ::Array{Polynomial{Int64,:y},2}) at .\REPL[40]:3
 [3] vcat(::Array{Polynomial{Int64,:x},2}, ::Array{Polynomial{Int64,:y},2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\SparseArrays\src\sparsevector.jl:1071
 [4] top-level scope at REPL[55]:1

julia> H = [[p p]; [q q]]
ERROR: ArgumentError: Polynomials have different indeterminates
Stacktrace:
 [1] assert_same_variable(::Symbol, ::Symbol) at C:\Users\Andreas\.julia\packages\Polynomials\1nhQ0\src\common.jl:394
 [2] promote_eltype(::Array{Polynomial{Int64,:x},2}, ::Array{Polynomial{Int64,:y},2}) at .\REPL[40]:3
 [3] vcat(::Array{Polynomial{Int64,:x},2}, ::Array{Polynomial{Int64,:y},2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\SparseArrays\src\sparsevector.jl:1071
 [4] top-level scope at REPL[56]:1

which are okay.

jverzani commented 3 years ago

One conversation on discourse caught my eye:

The new method:


function Base.promote_eltype(p::PP, q::QQ) where {T,X,N, P<:Polynomial{T,X}, PP<:Array{P,N},
                                                  S,Y,   Q<:Polynomial{S,Y}, QQ<:Array{Q,N}} 
    Polynomials.assert_same_variable(X,Y)
    Polynomial{promote_type(T,S),X}
end

Is essentially introducing something that is not otherwise expected for working with arrays arrays. (I also had to edit what I wrote above).

I'm not sure it matches what you want or not. (I had naively expected matrices constructed from sub-blocks to have promoted entries, but this is not the case unless a scalar is used for a 1x1 block.)

jverzani commented 3 years ago

Boy this is fiddly,

You found [p q 1] promotes, but [1 p q] does not, throwing an error. Seems like that needs to be corrected. Similarly with [p;q;1].

The last two (G = [[p p]; [one(p) one(q)]] for example) need this this extra method to throw that error. Unfortunately, G = [[one(q) one(p)]; [p p]] doesn't error (as the symbol for [one(q) one(p)] is :x

jverzani commented 3 years ago

Okay, getting some progress on this. Does the following work as expected on the patterns you would like to support? If so, I can generalize to other types and integrate into Polynomials.jl. Thanks for looking at this!

# how to integrate Polynomials.jl with matrices

using Polynomials
AbstractPolynomial = Polynomials.AbstractPolynomial

using Test
using LinearAlgebra

## These are modifications that need to be added to Polynomials.jl to get
## matrices to behave as expected.

_isconstant(p::Polynomials.AbstractPolynomial) = Polynomials.isconstant(p)
_isconstant(x) = true
_constantterm(p::Polynomials.AbstractPolynomial) = Polynomials.constantterm(p)
_constantterm(x) = x
_adjust_constant(p) = _isconstant(p) ? _constantterm(p) : p

function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:Polynomials.AbstractPolynomial{T}}
    Polynomials.isconstant(p) && return convert(S, Polynomials.constantterm(p))
    throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
end

Base.eltypeof(p::P) where {P <: Polynomials.AbstractPolynomial} = P
Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
                                               S,   Q<:AbstractPolynomial{S,X}} =
                                                   Polynomial{promote_type(T, S),X}

Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
                                              S,Y, Q<:AbstractPolynomial{S,Y}} =
                                                  Polynomials.assert_same_variable(X,Y)

function Base.promote_eltypeof(p::P, qs...) where {T,X, P<: Polynomials.AbstractPolynomial{T,X}}
    p′ = _adjust_constant(p)
    qs′ = _adjust_constant.(qs)
    Base.promote_type(Base.eltypeof(p′), Base.promote_eltypeof(qs′...))
end

function Base.promote_eltypeof(p::P, q::Q) where {T,X, P<: Polynomials.AbstractPolynomial{T,X},
                                                  S,Y, Q<: Polynomials.AbstractPolynomial{S,Y}}
    _isconstant(p) || _isconstant(q) || Polynomials.assert_same_variable(X,Y)
    p′ = _adjust_constant(p)
    q′ = _adjust_constant(q)
    Base.promote_type(Base.eltypeof(p′), Base.promote_eltypeof(q′))
end

# some testing
p = Polynomial([1,2],:x)
q = Polynomial([1,2],:y)
r = ImmutablePolynomial([1,2],:x)
s = ImmutablePolynomial([1,2],:y)

[p p]
[p r]
[p 1]
[1 1 p]
@test_throws ArgumentError [p q]
@test_throws ArgumentError [p s]
[p one(q)]
[p one(s)]
[p one(q) 1]
[1 one(q)]
[1 1 one(q)]
[1 one(q) p]
[1 1; one(q) p]
hcat([p p; 1 one(q)], I)

@test_throws ArgumentError [[p p]; [q q]]
[[p p]; [one(p) one(q)]]

andreasvarga commented 3 years ago

Sorry not being able to test earlier, but I was busy with the transition of DescriptorSystems to v2.0 of Polynomials. Finally, I succeded to update to a new release. The compilation times are once again longer (between 40 and 60 minutes), but this time I was not able to discover where things are jamming.

The above code is OK. But ...

I wonder if

julia> [p q I]
1×3 Array{Any,2}:
 Polynomial(1 + 2*x)  Polynomial(1 + 2*y)  UniformScaling{Bool}(true)

should produce an error, or the following error message is the right one

julia> [p q;I]
ERROR: ArgumentError: number of columns of each array must match (got (2, 1))
Stacktrace:
 [1] _typed_vcat(::Type{Any}, ::Tuple{Array{Any,2},Array{Any,2}}) at .\abstractarray.jl:1439
 [2] typed_vcat(::Type{Any}, ::Array{Any,2}, ::Array{Any,2}) at .\abstractarray.jl:1453
 [3] typed_hvcat(::Type{Any}, ::Tuple{Int64,Int64}, ::Polynomial{Int64,:x}, ::Vararg{Any,N} where N) at .\abstractarray.jl:1830
 [4] hvcat(::Tuple{Int64,Int64}, ::Polynomial{Int64,:x}, ::Vararg{Any,N} where N) at .\abstractarray.jl:1805
 [5] top-level scope at REPL[63]:1

Probably the following should work

julia> [p p;I]
ERROR: ArgumentError: number of columns of each array must match (got (2, 1))
Stacktrace:
 [1] _typed_vcat(::Type{Any}, ::Tuple{Array{Any,2},Array{Any,2}}) at .\abstractarray.jl:1439
 [2] typed_vcat(::Type{Any}, ::Array{Any,2}, ::Array{Any,2}) at .\abstractarray.jl:1453
 [3] typed_hvcat(::Type{Any}, ::Tuple{Int64,Int64}, ::Polynomial{Int64,:x}, ::Vararg{Any,N} where N) at .\abstractarray.jl:1830
 [4] hvcat(::Tuple{Int64,Int64}, ::Polynomial{Int64,:x}, ::Vararg{Any,N} where N) at .\abstractarray.jl:1805
 [5] top-level scope at REPL[66]:1

jverzani commented 3 years ago

Thanks!

The issue with [p q I] comes down to Base.promote_eltypeof(q,I) returning type Any, which matches what happens if q=1, say (e.g. [1 1.0 I] is similar, having eltype Any. I would guess that is how it should be.

Similarly, [p q;I] is an error through the same mechanism as [1 1.0; I]. That same error comes about with [p p; I]. Maybe such things could be worked around by adding methods to Base.promote_eltypeof for Arrays of Polynomials. Somewhere in this thread I did that.

But it seems that the right model for polynomials as matrix elements would be to treat non-constant polynomials as scalars of type P{T,X}; constant polynomials as scalars of type T, and not allow mixing of indeterminates under this view point. And that is all with the modifications to the underlying generic Julia code.

I did notice that what I had above needs even more! Basically [a,b] calls vect which calls Base.promote_typeof; whereas [a b] and [a;b] call cat which calls Base.promote_eltypeof. I'd guess this is related to how blocks are treated in matrix construction. At sometime you had sent in many modifications you had made to get this all to work and I was hopeful it could be done with much less work. Now I'm a bit afraid to go back and check :) For now I have similar definitions in each, of the following type:

_adjust_constant(p::Polynomials.AbstractPolynomial) = Polynomials.isconstant(p) ? Polynomials.constantterm(p) : p
_adjust_constant(x) = x

## [a,b] calls `Base.vect` which in turn calls Base.promote_typeof for promotion
Base.promote_typeof(p::P) where {P <: Polynomials.AbstractPolynomial} = Base.eltypeof(_adjust_constant(p))
function Base.promote_typeof(p::P, xs...) where {P <: Polynomials.AbstractPolynomial}
    x = _adjust_constant(p)
    Base.promote_type(Base.typeof(x), Base.promote_typeof(xs...))
end
function Base.promote_typeof(p::P, q::Q) where {T,X, P<: Polynomials.AbstractPolynomial{T,X},
                                                S,Y, Q<: Polynomials.AbstractPolynomial{S,Y}}
    isconstant(p) || isconstant(q) || Polynomials.assert_same_variable(X,Y)
    p′ = _adjust_constant(p)
    q′ = _adjust_constant(q)
    Base.promote_type(Base.eltypeof(p′), Base.eltypeof(q′))
end

## [p q] and [p;q] call `Base.cat` which in turn call Base.promote_eltypeof for promtion
Base.eltypeof(p::P) where {P <: Polynomials.AbstractPolynomial} = P
function Base.promote_eltypeof(p::P, xs...) where {T,X, P<: Polynomials.AbstractPolynomial{T,X}}
    x = _adjust_constant(p)
    Base.promote_type(Base.eltypeof(x), Base.promote_eltypeof(xs...))
end

function Base.promote_eltypeof(p::P, q::Q) where {T,X, P<: Polynomials.AbstractPolynomial{T,X},
                                                  S,Y, Q<: Polynomials.AbstractPolynomial{S,Y}}
    isconstant(p) || isconstant(q) || Polynomials.assert_same_variable(X,Y)
    p′ = _adjust_constant(p)
    q′ = _adjust_constant(q)
    Base.promote_type(Base.eltypeof(p′), Base.eltypeof(q′))
end

jverzani commented 3 years ago

I'm guessing, but just looked at DescriptorSystems. In RationalFunction.jl, the struct:

struct RationalTransferFunction{T} <: AbstractRationalTransferFunction
    num::Polynomial{T}        # numerator polynomial
    den::Polynomial{T}        # denominator polynomial
    ...

Should include the X, as in

struct RationalTransferFunction{T,X} <: AbstractRationalTransferFunction
    num::Polynomial{T,X}        # numerator polynomial
    den::Polynomial{T,X}
    ...

Unfortunately, this won't work with Polynomials v1.2; so if you were supporting both; you would want to define this conditionally on the version of Polynomials.

andreasvarga commented 3 years ago

Thanks. I will make the changes (I expect a small chain reaction). Still I wonder why the present version works, even with the "truncated" type. Do you have perhaps an answer?

john verzani @.***> schrieb am Fr., 19. März 2021, 16:24:

I'm guessing, but just looked at DescriptorSystems. In RationalFunction.jl, the struct:

struct RationalTransferFunction{T} <: AbstractRationalTransferFunction num::Polynomial{T} # numerator polynomial den::Polynomial{T} # denominator polynomial ...

Should include the X, as in

struct RationalTransferFunction{T,X} <: AbstractRationalTransferFunction num::Polynomial{T,X} # numerator polynomial den::Polynomial{T,X} ...

Unfortunately, this won't work with Polynomials v1.2; so if you were supporting both; you would want to define this conditionally on the version of Polynomials.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/312#issuecomment-802913167, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHECSBC2GADPJJ2ZOK5DTENUETANCNFSM4ZCNVYTQ .

jverzani commented 3 years ago

It isn't wrong, it just isn't as performant as having a fully realized type in the field of the struct. (If you don't add that X, you don't have a concrete type:

julia> Base.isconcretetype(Polynomial{Float64,X} where {X})
false

julia> Base.isconcretetype(Polynomial{Float64,:x})
true

) My guess is the slowdown occurs here. (But it is only a guess.)

jverzani commented 3 years ago

Sorry, want to pick your brain, as you have more developed ideas about what is expected than I do.

I have some promotion rules for mixing polynomials with constant polynomials and scalars that I think express this:

should promote to Polynomial type if mixed (except for Laurent)
non-constant polynomials must share the same indeterminate
constant polynomials of type P{T,X} should be treated as a constant of type T.
numeric constants should promote to a polynomial, when mixed

This is kinda working. Things like [1 p] promote to [Polynomial(1) p], for example when p is non-constant. But they have this side effect [1 one(p)] promotes to [1 1], which kind of makes sense, but I have a choice as to what becomes of [one(p)]. WIll it be [Polynomial(1)] or [1]? I could see either are being considered correct.

andreasvarga commented 3 years ago

Regarding the last issues, I obtained

julia> [1 one(Polynomial)]
1×2 Array{Polynomial{Float64,:x},2}:
 Polynomial(1.0)  Polynomial(1.0)

which is correct (in my opinion) (see also [1 one(im)]).

And also

julia> [one(Polynomial)]
1-element Array{Polynomial{Float64,:x},1}:
 Polynomial(1.0)

which is what I expected (compare with[one(1+im)]).

jverzani commented 3 years ago

In #328 I approach this problem a different way. I've been mumbling about what to do here for awhile, and haven't been really happy with what I had sketched out in #325 as it was too complicated and non-performant. In #328 a more explicit approach is required for mixing of variables. It may seem more cumbersome, but eliminates the question of what should [one(Polynomial,:x) one(Polynomial,:y)] be.

Would that be suitable for your needs?

andreasvarga commented 3 years ago

Sorry, I lost the point. I hope you will take the right decision.

john verzani @.***> schrieb am Mo., 29. März 2021, 21:49:

In #328 https://github.com/JuliaMath/Polynomials.jl/pull/328 I approach this problem a different way. I've been mumbling about what to do here for awhile, and haven't been really happy with what I had sketched out in #325 https://github.com/JuliaMath/Polynomials.jl/pull/325 as it was too complicated and non-performant. In #328 https://github.com/JuliaMath/Polynomials.jl/pull/328 a more explicit approach is required for mixing of variables. It may seem more cumbersome, but eliminates the question of what should [one(Polynomial,:x) one(Polynomial,:y)] be.

Would that be suitable for your needs?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/312#issuecomment-809665011, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEDAJ4DJZKF4CUNT6TTTGDKSXANCNFSM4ZCNVYTQ .

jverzani commented 3 years ago

Sorry, the basic point is for performance reasons the solution I want to merge requires the user to be specific about the element type when mixing terms with different variables through the use of a typed constructor.

So if p=Polynomial([1,2],:x) and q = Polynomial(3, :y) (a constant) we would have

[p,q] would throw an error (as :x is not :y)

Polynomial{Int,:x}[p,q] would be allowed, with q promoting to the specified type (as a constant polynomial can be converted to a polynomial with a different indeterminate).

Of course, in a function you would likely use typeof(p)[p, q] or some variant to get the type of p.

The use of numbers, which slowed your code down, would not need such annotations (e.g. [1,p] would just work).

andreasvarga commented 3 years ago

Thanks for the explanations. It's OK for me. I guess for rational functions the same approach will be used.

john verzani @.***> schrieb am Mo., 29. März 2021, 23:55:

Sorry, the basic point is for performance reasons the solution I want to merge requires the user to be specific about the element type when mixing terms with different variables through the use of a typed constructor.

So if p=Polynomial([1,2],:x) and q = Polynomial(3, :y) (a constant) we would have

[p,q] would throw an error (as :x is not :y)

Polynomial{Int,:x}[p,q] would be allowed, with q promoting to the specified type (as a constant polynomial can be converted to a polynomial with a different indeterminate).

Of course, in a function you would likely use typeof(p)[p, q] or some variant to get the type of p.

The use of numbers, which slowed your code down, would not need such annotations (e.g. [1,p] would just work).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/312#issuecomment-809739540, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEEZBUFQCT2OQIQN6HDTGDZM3ANCNFSM4ZCNVYTQ .

jverzani commented 3 years ago

Closed with #328

JuliaMath / Polynomials.jl

Mixing of variables in concatenation operations allowed #312