Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Just a question #308

Open
andreasvarga opened this issue Feb 5, 2021 · 30 comments
Open

Just a question #308

andreasvarga opened this issue Feb 5, 2021 · 30 comments

Comments

@andreasvarga
Copy link

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
Copy link
Author

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.

@jverzani
Copy link
Member

jverzani commented Feb 5, 2021

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.

@andreasvarga
Copy link
Author

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.

@jverzani
Copy link
Member

jverzani commented Feb 5, 2021

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).

@andreasvarga
Copy link
Author

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 andreasvarga mentioned this issue Feb 6, 2021
Open
@andreasvarga
Copy link
Author

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.

@jverzani
Copy link
Member

jverzani commented Feb 8, 2021

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.

@andreasvarga
Copy link
Author

andreasvarga commented Feb 9, 2021 via email

@andreasvarga
Copy link
Author

Another try.

polynomial_concatenations.zip

@andreasvarga
Copy link
Author

Please use this version
polynomial_concatenations.zip

@jverzani
Copy link
Member

jverzani commented Feb 9, 2021 via email

@jverzani
Copy link
Member

jverzani commented Feb 9, 2021 via email

@andreasvarga
Copy link
Author

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
Copy link
Author

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?

@jverzani
Copy link
Member

jverzani commented Feb 9, 2021 via email

@jverzani jverzani mentioned this issue Feb 9, 2021
Merged
10 tasks
@jverzani
Copy link
Member

jverzani commented Feb 9, 2021

#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.

@andreasvarga
Copy link
Author

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.

@jverzani
Copy link
Member

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 :)

@andreasvarga
Copy link
Author

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

@jverzani
Copy link
Member

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.

@andreasvarga
Copy link
Author

andreasvarga commented Feb 10, 2021 via email

@jverzani
Copy link
Member

jverzani commented Feb 11, 2021 via email

@andreasvarga
Copy link
Author

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
Copy link
Author

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?

@jverzani
Copy link
Member

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
Copy link
Member

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

@andreasvarga
Copy link
Author

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)

@jverzani
Copy link
Member

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)

@andreasvarga
Copy link
Author

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.

@jverzani
Copy link
Member

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 jverzani mentioned this issue Apr 12, 2021
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@jverzani @andreasvarga