Polynomials.jl icon indicating copy to clipboard operation
Polynomials.jl copied to clipboard

Just a question

Open andreasvarga opened this issue 4 years ago • 30 comments

I would appreciate very much if you could reveal how you enforce the polynomial type of elements in a polynomial matrix construct, as for example in the following vertical concatenation example:


julia> s = Polynomial([0,1],:s)
Polynomial(s)

julia> [s;0]
2-element Array{Polynomial{Int64},1}:
 Polynomial(s)
 Polynomial(0)

I am trying to implement a RationalTransferFunction type, where in a vertical concatenation I am getting an array of Any (instead an array of RationalTransferFunction), as below:

julia> s = rtf('s')
RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)

julia> [s;0]
2-element Array{Any,1}:
  RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)
 0

Many thank in advance for your advice.

andreasvarga avatar Feb 05 '21 13:02 andreasvarga

I have a mean to obtain this

rtf.([s;0])
2-element Array{RationalTransferFunction{Int64},1}:
 RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)
 RationalTransferFunction{Int64}(Polynomial(0), Polynomial(1), 0.0)

but I would prefer the elegant way you chose.

andreasvarga avatar Feb 05 '21 13:02 andreasvarga

Yeah it is a good question. I suspect that you get Any because the variable for s and the promotion of 0 to a polynomial will be different leading to an issue in the rtf function producing Any. Does the same thing happen if s uses the default variable :x? Is the rtf function in a package somewhere I could look at?

For the first concatenation [s;0] I would have to double check, but my guess is somewhere there is a call to promote_type(typeof(s), typeof(0)) which would return Polynomial. However, since the variable is a field and not in the type, it will take the default variable name :x. That's more than a bit of a nuisance. Makes me think a redesign of the basic type might be called for to solve this.

jverzani avatar Feb 05 '21 13:02 jverzani

The RationalTransferFunction will be a new object in the next release of the DescriptorSystems.jl. My hope was that somebody will implement a rational function package, which I can use as basis for my transfer function object. I remember that I already mentioned this issue earlier to you. However, apparently the existing RationalFunctions is not further developed, so I decided to implement this object and related functions inspired by several other packages, including Polynomials. My hope was to discover there a function for concatenation of elements, but I was not able to find one. So, I dared to put you my question why

julia> [s 1.5; 1 1]
2×2 Array{Polynomial{Float64},2}:
 Polynomial(1.0*s)  Polynomial(1.5)
 Polynomial(1.0)    Polynomial(1.0)

works. Certainly, I can try to implement a complete set of concatenation function (i.e., hcat, vcat, hvcat, etc.), but this exercise was already very demanding for concatenating systems with matrices and UniformScaling, so my hope was to learn from your tricks. I look forward if you have an explanation for the above behaviour.

And, of course, thanks for the fast reply.

andreasvarga avatar Feb 05 '21 14:02 andreasvarga

Yes, concatenation with Polynomial objects falls back to Julia's default, which means you have to get the promotion mechanisms set up as you want them. (This is why my suspicion about variables being different is some cause of your "Any" type). The promotion is handled in the abstract.jl file (https://github.com/JuliaMath/Polynomials.jl/blob/master/src/abstract.jl#L44).

jverzani avatar Feb 05 '21 15:02 jverzani

Thanks for the hint. I implemented the promotion rule, but it still remain two issues. The first one concerns the possibility of different variables (this is an issue for Polynomials too):

julia> s = Polynomial([0,1],:s)
Polynomial(s)

julia> z = Polynomial([0,1],:z)
Polynomial(z)

julia> [s;z]
2-element Array{Polynomial{Int64},1}:
 Polynomial(s)
 Polynomial(z)

The second concerns the possibility of different sampling time (0 stands for continuous-time, positive or -1 stand for discrete-time): 

julia> s = rtf('s')
RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)

julia> z = rtf('z',Ts = 1)
RationalTransferFunction{Int64}(Polynomial(z), Polynomial(1), 1.0)

julia> [s;z]
2-element Array{RationalTransferFunction{Int64},1}:
 RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)
 RationalTransferFunction{Int64}(Polynomial(z), Polynomial(1), 1.0)

I set the sampling time -Inf for undefined sampling time, in the hope to be able to handle the above issue later, so that now:

julia> [s;0]
2-element Array{RationalTransferFunction{Int64},1}:
 RationalTransferFunction{Int64}(Polynomial(s), Polynomial(1), 0.0)
 RationalTransferFunction{Int64}(Polynomial(0), Polynomial(1), -Inf)

Unfortunately, I see no simple solution. I must probably implement a new set of concatenation functions.

andreasvarga avatar Feb 05 '21 17:02 andreasvarga

I implemented a set of concatenation functions for polynomial matrices. It is now possible to generate polynomial matrices with correctly set variable names. For example, the following examples combines scalars, matrices and UniformScalings into one matrix, with the type automatically adjusted and the variable names correctly set:

julia> s = Polynomial([0, 1],'s'); z = Polynomial([0, 1],'z');

julia> Gc = [s^2 1; -s 2; I; rand(1,2)]
5×2 Array{Polynomial{Float64},2}:
 Polynomial(1.0*s^2)   Polynomial(1.0)
 Polynomial(-1.0*s)    Polynomial(2.0)
 Polynomial(1.0)       Polynomial(0.0)
 Polynomial(0.0)       Polynomial(1.0)
 Polynomial(0.127124)  Polynomial(0.40111)

julia> variable.(Gc)
5×2 Array{Polynomial{Float64},2}:    
 Polynomial(1.0*s)  Polynomial(1.0*s)
 Polynomial(1.0*s)  Polynomial(1.0*s)
 Polynomial(1.0*s)  Polynomial(1.0*s)
 Polynomial(1.0*s)  Polynomial(1.0*s)
 Polynomial(1.0*s)  Polynomial(1.0*s)

julia> [s;z]
ERROR: all polynomial matrix elements must have the same variable
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] promote_poly_var(::Polynomial{Int64}, ::Vararg{Polynomial{Int64},N} where N) at C:\Users\Andreas\Documents\software\Julia\DescriptorSystems.jl\src\polynomial_concatenations.jl:235
 [3] vcat(::Polynomial{Int64}, ::Polynomial{Int64}) at C:\Users\Andreas\Documents\software\Julia\DescriptorSystems.jl\src\polynomial_concatenations.jl:22
 [4] top-level scope at REPL[9]:1

If you want to include this stuff in the Polynomials as additional support to work with matrices, I would be pleased to forward you the software.

andreasvarga avatar Feb 08 '21 17:02 andreasvarga

Neat. I'm hoping to make this easier for you. I'm almost done with a PR that moves the polynomial symbol into the type parameter so things like promote will just work. It took a fair amount of code churn, so might be a bit of time before it is tagged, but I should have something to play around with by tomorrow. I'd love to see what you put together for this task on your end, as hopefully I can either borrow or streamline what is necessary.

jverzani avatar Feb 08 '21 21:02 jverzani

I hope it works.

john verzani [email protected] schrieb am Mo., 8. Feb. 2021, 22:31:

Neat. I'm hoping to make this easier for you. I'm almost done with a PR that moves the polynomial symbol into the type parameter so things like promote will just work. It took a fair amount of code churn, so might be a bit of time before it is tagged, but I should have something to play around with by tomorrow. I'd love to see what you put together for this task on your end, as hopefully I can either borrow or streamline what is necessary.

— 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/308#issuecomment-775476307, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEHSILEG6HFKARK5URLS6BJ2JANCNFSM4XEY24RQ .

andreasvarga avatar Feb 09 '21 09:02 andreasvarga

Another try.

polynomial_concatenations.zip

andreasvarga avatar Feb 09 '21 09:02 andreasvarga

Please use this version polynomial_concatenations.zip

andreasvarga avatar Feb 09 '21 12:02 andreasvarga

Thanks. Before I dive in, do you have some test suite?

On Tue, Feb 9, 2021 at 7:15 AM Andreas Varga [email protected] wrote:

Please use this version polynomial_concatenations.zip https://github.com/JuliaMath/Polynomials.jl/files/5951254/polynomial_concatenations.zip

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/308#issuecomment-775895275, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TGXWHHDZE3MFOSJR5LS6ERO5ANCNFSM4XEY24RQ .

-- John Verzani Department of Mathematics College of Staten Island, CUNY [email protected]

jverzani avatar Feb 09 '21 12:02 jverzani

Just to be clear, with my current branch, I can do all this before adding your modifications, but am not sure what more you are seeking:


julia> p = Polynomial([1,2,3], :y)

Polynomial(1 + 2*y + 3*y^2)


julia> A = [1 p; p^2 p]

2×2 Array{Polynomial{Int64,:y},2}:

 Polynomial(1)                                  Polynomial(1 + 2*y + 3*y^2)

 Polynomial(1 + 4*y + 10*y^2 + 12*y^3 + 9*y^4)  Polynomial(1 + 2*y + 3*y^2)


julia> vcat(A,I)

4×2 Array{Polynomial{Int64,:y},2}:

 Polynomial(1)                                  Polynomial(1 + 2*y + 3*y^2)

 Polynomial(1 + 4*y + 10*y^2 + 12*y^3 + 9*y^4)  Polynomial(1 + 2*y + 3*y^2)

 Polynomial(1)                                  Polynomial(0)

 Polynomial(0)                                  Polynomial(1)


julia> hcat(A,I)

2×4 Array{Polynomial{Int64,:y},2}:

 Polynomial(1)                                  …  Polynomial(1)
Polynomial(0)

 Polynomial(1 + 4*y + 10*y^2 + 12*y^3 + 9*y^4)     Polynomial(0)
Polynomial(1)


julia> hvcat((2,2), A,I,A,I)

4×4 Array{Polynomial{Int64,:y},2}:

 Polynomial(1)                                  …  Polynomial(1)
Polynomial(0)

 Polynomial(1 + 4*y + 10*y^2 + 12*y^3 + 9*y^4)     Polynomial(0)
Polynomial(1)

 Polynomial(1)                                     Polynomial(1)
Polynomial(0)

 Polynomial(1 + 4*y + 10*y^2 + 12*y^3 + 9*y^4)     Polynomial(0)
Polynomial(1)

But this is failing:


julia> A + I

ERROR: MethodError: convert(::Type{Union{}}, ::Polynomial{Int64,:y}) is
ambiguous.

On Tue, Feb 9, 2021 at 7:36 AM John Verzani [email protected] wrote:

Thanks. Before I dive in, do you have some test suite?

On Tue, Feb 9, 2021 at 7:15 AM Andreas Varga [email protected] wrote:

Please use this version polynomial_concatenations.zip https://github.com/JuliaMath/Polynomials.jl/files/5951254/polynomial_concatenations.zip

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/308#issuecomment-775895275, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TGXWHHDZE3MFOSJR5LS6ERO5ANCNFSM4XEY24RQ .

-- John Verzani Department of Mathematics College of Staten Island, CUNY [email protected]

-- John Verzani Department of Mathematics College of Staten Island, CUNY [email protected]

jverzani avatar Feb 09 '21 12:02 jverzani

Apparently you can do everything what I can also do! Or perhaps not ?

try [ A rand(2,2);I;[1 2p im 3.]]

Just in case you would like to have a look at the last version: polynomial_concatenations.zip

andreasvarga avatar Feb 09 '21 14:02 andreasvarga

Unfortunately, I have now the problem:

[zeros(2,2) I]
2×4 Array{Polynomial{Float64},2}:
 Polynomial(0.0)  Polynomial(0.0)  Polynomial(1.0)  Polynomial(0.0)
 Polynomial(0.0)  Polynomial(0.0)  Polynomial(0.0)  Polynomial(1.0)

Do you know a trick to avoid executing my code and falling back to the standard linear algebra code?

andreasvarga avatar Feb 09 '21 14:02 andreasvarga

Just wait :)

I should have these changes on a branch soon that you can explore. Since it is breaking (technically) I want to get in a few other breaking changes before I tag.

As for your last example, I still have some work to do:


julia> [ A rand(2,2);I;[1 2p im 3.]]

ERROR: TypeError: in typeassert, expected Polynomial{Complex{Float64},:y},
got a value of type Polynomial{Int64,:y}

On Tue, Feb 9, 2021 at 9:42 AM Andreas Varga [email protected] wrote:

Unfortunately, I have now the problem:

[zeros(2,2) I]

2×4 Array{Polynomial{Float64},2}:

Polynomial(0.0) Polynomial(0.0) Polynomial(1.0) Polynomial(0.0)

Polynomial(0.0) Polynomial(0.0) Polynomial(0.0) Polynomial(1.0)

Do you know a trick to avoid executing my code and falling back to the standard linear algebra code?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/308#issuecomment-775990202, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TAT7TZLBFCNJBUWFV3S6FCWXANCNFSM4XEY24RQ .

-- John Verzani Department of Mathematics College of Staten Island, CUNY [email protected]

jverzani avatar Feb 09 '21 14:02 jverzani

#310 contains a PR that moves the symbol into a type parameter. It should help clean up what you need quite a bit (and maybe reduce your modifications to none). Let me know how it goes.

jverzani avatar Feb 09 '21 19:02 jverzani

I finished with my implementation of Polynomial concatenations. This was a good exercise which now serves as basis of implementing similar software for rational transfer functions. Such an object I defined as:

abstract type AbstractRationalTransferFunction end
struct RationalTransferFunction{T} <: AbstractRationalTransferFunction
    num::Polynomial{T}          # numerator polynomial
    den::Polynomial{T}           # denominator polynomial
    Ts::Union{Real,Nothing}   # sampling time
    function RationalTransferFunction{T}(num::Polynomial{T}, den::Polynomial{T}, Ts::Union{Real,Nothing}) where T <: Number
        length(num) > 1 && length(den) > 1 && num.var != den.var && 
              error("The numerator and denominator polynomials must have the same variable")
        if all(den == zero(den))
            error("Cannot create a rational function with zero denominator")
        elseif all(num == zero(num))
            # The numerator is zero, make the denominator 1
            den = one(den)
        end
        # Validate sampling time
        isnothing(Ts) || Ts >= 0 || Ts == -1 || 
             error("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")
        new{T}(num, den, Ts)
    end
end

The sampling time Ts is a system specific quantity and Ts = nothing would correspond to pure rational functions. I wonder if the parameter Ts can also be moved to the object type, as you did with the polynom variable.

Could you tell me what impact will have your changes on my (an other) software relying on Polynomials. Perhaps none?

PS. I don't know how to use #310. Probably, the best is for me to wait (patiently) the final result of your modification.

andreasvarga avatar Feb 10 '21 13:02 andreasvarga

Thanks. I don't have time right now for review, but on first glance this is how I would approach it. (One thing the construct den.var will not work with the new version, rather Polynomials.indeterminate(den) is the new idiom.

As for moving Ts into the type, I'd think that would be an issue if there were many possible values, as it would mean different types for each problem, and hence many recompilations. For the symbol that isn't really the expected case. Though if there were a few, then sure that might be useful for the same reason moving the symbol up to that level is.

As for using #310, I think you can:

  • create a new directory
  • activate that directory: ] activate .
  • add the branch: ] add https://github.com/jverzani/Polynomials.jl#v2.0.0
  • then the usual using Polynomials

(hopefully :)

jverzani avatar Feb 10 '21 15:02 jverzani

It works ! Congratulations!

I wonder if you must use _indeterminate(PP) === nothing instead _indeterminate(PP) == nothing below (in common.jl)

function indeterminate(PP::Type{P}, p::AbstractPolynomial) where {P <: AbstractPolynomial}
    X = _indeterminate(PP) == nothing ? indeterminate(p) :  _indeterminate(PP)
end

andreasvarga avatar Feb 10 '21 16:02 andreasvarga

I might have this wrong, but I think == is enough, thought I have seen === used as well. I can certainly change it, I don't see the harm and perhaps there are benefits I'm unaware of.

Let me know if I got this setup correctly for your use. It is a lot of code churn and hopefully makes your task much simpler if not ready made.

jverzani avatar Feb 10 '21 19:02 jverzani

Would you suggest to modify my codes to use the new Polynomials? When do you think the changes go into a next release ?

john verzani [email protected] schrieb am Mi., 10. Feb. 2021, 20:20:

I might have this wrong, but I think == is enough, thought I have seen === used as well. I can certainly change it, I don't see the harm and perhaps there are benefits I'm unaware of.

Let me know if I got this setup correctly for your use. It is a lot of code churn and hopefully makes your task much simpler if not ready made.

— 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/308#issuecomment-776954270, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEBMNDSZC34JCDOWVVTS6LL7JANCNFSM4XEY24RQ .

andreasvarga avatar Feb 10 '21 20:02 andreasvarga

Hopefully in a week or two I’ll tag a release.

On Wed, Feb 10, 2021 at 3:59 PM Andreas Varga [email protected] wrote:

Would you suggest to modify my codes to use the new Polynomials? When do you think the changes go into a next release ?

john verzani [email protected] schrieb am Mi., 10. Feb. 2021, 20:20:

I might have this wrong, but I think == is enough, thought I have seen === used as well. I can certainly change it, I don't see the harm and perhaps there are benefits I'm unaware of.

Let me know if I got this setup correctly for your use. It is a lot of code churn and hopefully makes your task much simpler if not ready made.

— 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/308#issuecomment-776954270 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/ALJDHEBMNDSZC34JCDOWVVTS6LL7JANCNFSM4XEY24RQ

.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/308#issuecomment-777031712, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TFHJDE74ILQEBK6EITS6LXS3ANCNFSM4XEY24RQ .

-- John Verzani Department of Mathematics College of Staten Island, CUNY [email protected]

jverzani avatar Feb 11 '21 13:02 jverzani

This is a short account of my first experince moving to the new Polynomials.

I started to modify the software in MatrixPencils.jl to comply with the new setup of Polynomials. I performed modifications trying to cover both the current version as well as the new one.

In some cases I just added a new call as below:

function pmdeg(P::Union{AbstractVecOrMat{Polynomial{T}},Polynomial{T},Number}) where T
   typeof(P) <: Number && (P == 0 ? (return -1) : (return 0) )
   return maximum(degree.(P))
end

function pmdeg(P::Union{AbstractVecOrMat{Polynomial{T,X}},Polynomial{T,X},Number}) where {T,X}
   typeof(P) <: Number && (P == 0 ? (return -1) : (return 0) )
   return maximum(degree.(P))
end

In other cases this duplication was not possible, since I received errors complainig of ambiguous methods, but the code works without any change! Apparently the type Matrix{Polynomial{T}} automatically covers Matrix{Polynomial{T,X}} .

And some cases I was not able to manage with either of the two approaches. For example, I got a polynomial matrix displayed as:


julia> Q
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(1.0)  Polynomial(0)
 Polynomial(0)    Polynomial(2.0)
julia> Q
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(1.0)  Polynomial(0)
 Polynomial(0)    Polynomial(2.0)

julia> Polynomials.indeterminate.(Q)
2×2 Array{Symbol,2}:
 :s  :s
 :s  :s

julia> eltype(Q)
Polynomial{T,:s} where T<:Number

julia> typeof(Q)
Array{Polynomial{T,:s} where T<:Number,2}

Question: I wonder if Q is a valid polynomial matrix, since I expected Q of type 2×2 Array{Polynomial{Float64,:s},2}.

So, I added a method for this type too, as below:

function poly2pm(PM::AbstractMatrix{Polynomial{T,X} where T <: Number}; grade::Union{Int,Missing} = missing) where X

The value of T I have to figure out using eltype(eltype(PM)). And, it works!

Is any systematics which you could recommed for complying with both old and new versions? In this moment, it seems to me to be a matter of trials and errors.

andreasvarga avatar Feb 11 '21 15:02 andreasvarga

With the new structure, we could consider the following structures:


julia> tx=AbstractVecOrMat{Polynomial{T,X} where T} where {X}
Union{AbstractArray{Polynomial{T,X} where T,1}, AbstractArray{Polynomial{T,X} where T,2}} where X

julia> tt =AbstractVecOrMat{Polynomial{T,X} where X} where {T}
Union{AbstractArray{Polynomial{T,X} where X,1}, AbstractArray{Polynomial{T,X} where X,2}} where T

julia> ttx = AbstractVecOrMat{Polynomial{T,X}} where {T,X}
Union{AbstractArray{Polynomial{T,X},1}, AbstractArray{Polynomial{T,X},2}} where X where T

julia> t = AbstractVecOrMat{Polynomial}
Union{AbstractArray{Polynomial,1}, AbstractArray{Polynomial,2}}

The first three are clearly different, but the last two are the same, because

julia> t <: ttx && ttx <: t
true

The question is: in the most general case of a polynomial matrix input the Union{tt,tx,ttx} should be always supported?

andreasvarga avatar Feb 11 '21 17:02 andreasvarga

This is a short account of my first experince moving to the new Polynomials.

I started to modify the software in MatrixPencils.jl to comply with the new setup of Polynomials. I performed modifications trying to cover both the current version as well as the new one.

In some cases I just added a new call as below:

function pmdeg(P::Union{AbstractVecOrMat{Polynomial{T}},Polynomial{T},Number}) where T
   typeof(P) <: Number && (P == 0 ? (return -1) : (return 0) )
   return maximum(degree.(P))
end

function pmdeg(P::Union{AbstractVecOrMat{Polynomial{T,X}},Polynomial{T,X},Number}) where {T,X}
   typeof(P) <: Number && (P == 0 ? (return -1) : (return 0) )
   return maximum(degree.(P))
end

In other cases this duplication was not possible, since I received errors complainig of ambiguous methods, but the code works without any change! Apparently the type Matrix{Polynomial{T}} automatically covers Matrix{Polynomial{T,X}} .

And some cases I was not able to manage with either of the two approaches. For example, I got a polynomial matrix displayed as:


julia> Q
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(1.0)  Polynomial(0)
 Polynomial(0)    Polynomial(2.0)
julia> Q
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(1.0)  Polynomial(0)
 Polynomial(0)    Polynomial(2.0)

julia> Polynomials.indeterminate.(Q)
2×2 Array{Symbol,2}:
 :s  :s
 :s  :s

julia> eltype(Q)
Polynomial{T,:s} where T<:Number

julia> typeof(Q)
Array{Polynomial{T,:s} where T<:Number,2}

Question: I wonder if Q is a valid polynomial matrix, since I expected Q of type 2×2 Array{Polynomial{Float64,:s},2}.

So, I added a method for this type too, as below:

function poly2pm(PM::AbstractMatrix{Polynomial{T,X} where T <: Number}; grade::Union{Int,Missing} = missing) where X

The value of T I have to figure out using eltype(eltype(PM)). And, it works!

Is any systematics which you could recommed for complying with both old and new versions? In this moment, it seems to me to be a matter of trials and errors.

How is Q constructed? The promotion mechanism isn't kicking in as expected and one polynomial has type {Int, :s}, the other type {Float64, :s} hence the unexpected type of the matrix.

jverzani avatar Feb 12 '21 01:02 jverzani

With the new structure, we could consider the following structures:


julia> tx=AbstractVecOrMat{Polynomial{T,X} where T} where {X}
Union{AbstractArray{Polynomial{T,X} where T,1}, AbstractArray{Polynomial{T,X} where T,2}} where X

julia> tt =AbstractVecOrMat{Polynomial{T,X} where X} where {T}
Union{AbstractArray{Polynomial{T,X} where X,1}, AbstractArray{Polynomial{T,X} where X,2}} where T

julia> ttx = AbstractVecOrMat{Polynomial{T,X}} where {T,X}
Union{AbstractArray{Polynomial{T,X},1}, AbstractArray{Polynomial{T,X},2}} where X where T

julia> t = AbstractVecOrMat{Polynomial}
Union{AbstractArray{Polynomial,1}, AbstractArray{Polynomial,2}}

The first three are clearly different, but the last two are the same, because

julia> t <: ttx && ttx <: t
true

The question is: in the most general case of a polynomial matrix input the Union{tt,tx,ttx} should be always supported?

Maybe you want t = AbstractVecOrMat{<: Polynomial} which should allow for the case of unspecified T, X

jverzani avatar Feb 12 '21 01:02 jverzani

The following commands lead to X, Q and R which escape the promotion mechanism:

s = Polynomial([0, 1],:s);
NUM = [s^2+3*s+3 1; -1 2*s^2+7*s+4];
DEN = [(s+1)^2 s+2; (s+1)^3 (s+1)*(s+2)];
X = divrem.(NUM,DEN)
Q = first.(X)
R = last.(X)

julia> X = divrem.(NUM,DEN)
2×2 Array{Tuple{Polynomial{T,:s} where T<:Number,Polynomial{T,:s} where T<:Number},2}:
 (Polynomial(1.0), Polynomial(2.0 + 1.0*s))  (Polynomial(0), Polynomial(1))
 (Polynomial(0), Polynomial(-1))             (Polynomial(2.0), Polynomial(1.0*s))

julia> Q = first.(X)
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(1.0)  Polynomial(0)
 Polynomial(0)    Polynomial(2.0)

julia> R = last.(X)
2×2 Array{Polynomial{T,:s} where T<:Number,2}:
 Polynomial(2.0 + 1.0*s)  Polynomial(1)
 Polynomial(-1)           Polynomial(1.0*s)

andreasvarga avatar Feb 12 '21 09:02 andreasvarga

Yes indeed. The issue is promote between containers with polynomials. In this example there is no promotion between Tuple{P{Int,:X},2) and Tuple{P{Float64,:X},2) which is understandable, tuples and all. But I naively expected that if I made those Vectors (in drem below), there would be promotion, but that is not the case and it takes a separate step.

I'm not sure where to do this, but this pattern works:

drem(u,v) = [divrem(u,v)...] # use vectors, not tuples
X = drem.(NUM,DEN)          # 2×2 Array{Array{T,1} where T,2}
convert(Array{eltype(promote(X...)),2},  X)

jverzani avatar Feb 12 '21 12:02 jverzani

Sorry disturbing you once again.

I started to update my current projects to use v2.0 of Polynomials. Great work! Congratulations!

During transition of MatrixPencils.jl, I got a message from the CompatHelper to change in Project.toml

[compat]
Polynomials = "2.0" 

to

[compat]
Polynomials = "1.2, 2.0"

in order to support previous versions working with Polynomials v1.2.

I wonder if this has any sense for me. I would appreciate learning what is your opinion in this respect and if this is also an issue for you, i.e., to comply with the previous versions inside Polynomials. Thank you for your time.

andreasvarga avatar Mar 12 '21 11:03 andreasvarga

I'll have to look. My guess would be that if you do one of three things some workaround will be needed:

  • subtype AbstractPolynomial{T} (in which case AbstractPolynomial{T,X} is needed in 2.0.0, but not < 2.0.0
  • use Polynomial{X} or some other type in a struct (in which case, it will be more performant to use Polynomial{T,X} in the struct in 2.0 but not < 2.0.0
  • access the indeterminate with p.var (in which case defining indeterminate(p)=p.var if the VERSION is < 2.0.0 would be needed)

I should be able to have a look next week through your package if the above isn't enough guidance.

jverzani avatar Mar 12 '21 16:03 jverzani