Updates for Polynomials v1.0 (#292)

* Updates for Polynomials v0.7

* fix test failures

Project.toml

@@ -27,5 +27,5 @@ IterTools = "1.0"
 LaTeXStrings = "1.0"
 OrdinaryDiffEq = "5.2"
 Plots = "0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 1.0"
-Polynomials = "0.6.0"
+Polynomials = "0.7, 0.8, 1.0"
 julia = "1.0"

example/symbolic_computations.jl

@@ -12,7 +12,7 @@ using SymPy
 Base.complex(::Type{SymPy.Sym}) = Sym
 Base.eigvals(a::Matrix{Sym}) = Vector{Sym}(collect(keys(call_matrix_meth(a, :eigenvals))))
 
-function Polynomials.roots(p::Polynomials.Poly{SymPy.Sym})
+function Polynomials.roots(p::Polynomials.Polynomial{SymPy.Sym})
     x = Sym("x")
     return SymPy.solve(p(x), x)
 end

src/ControlSystems.jl

@@ -90,7 +90,7 @@ export  LTISystem,
 # QUESTION: are these used? LaTeXStrings, Requires, IterTools
 using Plots, LaTeXStrings, LinearAlgebra
 import Polynomials
-import Polynomials: Poly, coeffs, polyval
+import Polynomials: Polynomial, coeffs
 using OrdinaryDiffEq, DelayDiffEq
 export Plots
 import Base: +, -, *, /, (==), (!=), isapprox, convert, promote_op

src/delay_systems.jl

@@ -87,13 +87,13 @@ end
 function _lsim(sys::DelayLtiSystem{T,S}, Base.@nospecialize(u!), t::AbstractArray{<:Real}, x0::Vector{T}, alg; kwargs...) where {T,S}
     P = sys.P
 
-    if ~iszero(P.D22)
+    if !iszero(P.D22)
         throw(ArgumentError("non-zero D22-matrix block is not supported")) # Due to limitations in differential equations
     end
 
     t0 = first(t)
     dt = t[2] - t[1]
-    if ~all(diff(t) .≈ dt) # QUESTION Does this work or are there precision problems?
+    if !all(diff(t) .≈ dt) # QUESTION Does this work or are there precision problems?
         error("The t-vector should be uniformly spaced, t[2] - t[1] = $dt.") # Perhaps dedicated function for checking this?
     end
 
@@ -219,25 +219,25 @@ end
 
 # Used for pade approximation
 """
-`p2 = _linscale(p::Poly, a)`
+`p2 = _linscale(p::Polynomial, a)`
 
 Given a polynomial `p` and a number `a, returns the polynomials `p2` such that
 `p2(s) == p(a*s)`.
 """
-function _linscale(p::Poly, a)
+function _linscale(p::Polynomial, a)
     # This function should perhaps be implemented in Polynomials.jl
-    coeffs_scaled = similar(p.a, promote_type(eltype(p), typeof(a)))
+    coeffs_scaled = similar(p.coeffs, promote_type(eltype(p), typeof(a)))
     a_pow = 1
-    coeffs_scaled[1] = p.a[1]
-    for k=2:length(p.a)
+    coeffs_scaled[1] = p.coeffs[1]
+    for k=2:length(p.coeffs)
         a_pow *= a
-        coeffs_scaled[k] = p.a[k]*a_pow
+        coeffs_scaled[k] = p.coeffs[k]*a_pow
     end
-    return Poly(coeffs_scaled)
+    return Polynomial(coeffs_scaled)
 end
 
 # Coefficeints for Padé approximations
-# PADE_Q_COEFFS = [Poly([binomial(N,i)*prod(N+1:2*N-i) for i=0:N]) for N=1:10]
+# PADE_Q_COEFFS = [Polynomial([binomial(N,i)*prod(N+1:2*N-i) for i=0:N]) for N=1:10]
 const PADE_Q_COEFFS =  [[2, 1],
  [12, 6, 1],
  [120, 60, 12, 1],
@@ -257,7 +257,7 @@ Compute the `N`th order Padé approximation of a time-delay of length `τ`.
 function pade(τ::Real, N::Int)
     if !(1 <= N <= 10); error("Order of Padé approximation must be between 1 and 10. Got $N."); end
 
-    Q = Poly(PADE_Q_COEFFS[N])
+    Q = Polynomials.Polynomial(PADE_Q_COEFFS[N])
 
     return tf(_linscale(Q, -τ), _linscale(Q, τ)) # return Q(-τs)/Q(τs)
 end

src/discrete.jl

@@ -166,15 +166,14 @@ function toeplitz(c,r)
     A
 end
 
-import Polynomials
 
 """
 `c2d_roots2poly(ro,h)`
 
 returns the polynomial coefficients in discrete time given a vector of roots in continuous time
 """
 function c2d_roots2poly(ro,h)
-    return real((Polynomials.poly(exp(ro.*h))).a[end:-1:1])
+    return real((Polynomials.poly(exp(ro.*h))).coeffs[end:-1:1])
 end
 
 """
@@ -183,8 +182,8 @@ end
 returns the polynomial coefficients in discrete time given polynomial coefficients in continuous time
 """
 function c2d_poly2poly(p,h)
-    ro = Polynomials.roots(Polynomials.Poly(p[end:-1:1]))
-    return real(Polynomials.poly(exp(ro.*h)).a[end:-1:1])
+    ro = Polynomials.roots(Polynomials.Polynomial(p[end:-1:1]))
+    return real(Polynomials.poly(exp(ro.*h)).coeffs[end:-1:1])
 end
 
 

src/synthesis.jl

@@ -132,7 +132,7 @@ end
 function acker(A,B,P)
     n = length(P)
     #Calculate characteristic polynomial
-    poly = mapreduce(p -> Poly([1, -p]), *, P, init=Poly(one(eltype(P))))
+    poly = mapreduce(p -> Polynomial([1, -p]), *, P, init=Polynomial(one(eltype(P))))
     q = zero(Array{promote_type(eltype(A),Float64),2}(undef, n,n))
     for i = n:-1:0
         q += A^(n-i)*poly[i]

src/types/SisoTfTypes/SisoRational.jl

@@ -1,8 +1,8 @@
 ## User should just use TransferFunction
 struct SisoRational{T} <: SisoTf{T}
-    num::Poly{T}
-    den::Poly{T}
-    function SisoRational{T}(num::Poly{T}, den::Poly{T}) where T <: Number
+    num::Polynomial{T}
+    den::Polynomial{T}
+    function SisoRational{T}(num::Polynomial{T}, den::Polynomial{T}) where T <: Number
         if all(den == zero(den))
             error("Cannot create SisoRational with zero denominator")
         elseif all(num == zero(num))
@@ -12,14 +12,14 @@ struct SisoRational{T} <: SisoTf{T}
         new{T}(num, den)
     end
 end
-function SisoRational(num::Poly{T1}, den::Poly{T2}) where T1 <: Number where T2 <: Number
+function SisoRational(num::Polynomial{T1}, den::Polynomial{T2}) where T1 <: Number where T2 <: Number
     T = promote_type(T1,T2)
-    SisoRational{T}(Poly{T}(num.a), Poly{T}(den.a))
+    SisoRational{T}(Polynomial{T}(num.coeffs), Polynomial{T}(den.coeffs))
 end
-SisoRational{T}(num::Poly, den::Poly) where T = SisoRational{T}(convert(Poly{T}, num), convert(Poly{T}, den))
+SisoRational{T}(num::Polynomial, den::Polynomial) where T = SisoRational{T}(convert(Polynomial{T}, num), convert(Polynomial{T}, den))
 
 function SisoRational{T}(num::AbstractVector, den::AbstractVector) where T <: Number # NOTE: Typearguemnts on the parameters?
-    SisoRational{T}(Poly{T}(reverse(num)), Poly{T}(reverse(den)))
+    SisoRational{T}(Polynomial{T}(reverse(num)), Polynomial{T}(reverse(den)))
 end
 function SisoRational(num::AbstractVector{T1}, den::AbstractVector{T2}) where T1 <: Number where T2 <: Number
     T = promote_type(T1,T2)
@@ -81,11 +81,11 @@ pole(f::SisoRational) = roots(f.den)
 
 function evalfr(f::SisoRational{T}, s::Number) where T
     S = promote_op(/, promote_type(T, typeof(s)), promote_type(T, typeof(s)))
-    den = polyval(f.den, s)
+    den = f.den(s)
     if den == zero(S)
         convert(S,Inf)
     else
-        polyval(f.num, s)/den
+        f.num(s)/den
     end
 end
 

src/types/SisoTfTypes/SisoZpk.jl

@@ -149,9 +149,9 @@ function evalfr(f::SisoZpk{T1,TR}, s::T2) where {T1<:Number, TR<:Number, T2<:Num
 end
 
 
-function poly_factors2string(poly_factors::AbstractArray{<:Poly{T}}, var) where T
+function poly_factors2string(poly_factors::AbstractArray{<:Polynomial{T}}, var) where T
     if length(poly_factors) == 0
-        str = sprint(printpolyfun(var), Poly(one(T)))
+        str = sprint(printpolyfun(var), Polynomial(one(T)))
     elseif length(poly_factors) == 1
         str = sprint(printpolyfun(var), poly_factors[1])
     else
@@ -233,7 +233,7 @@ function +(f1::SisoZpk{T1,TR1}, f2::SisoZpk{T2,TR2}) where {T1<:Number,T2<:Numbe
     # FIXME:
     # Threshold for pole-zero cancellation should depend on the roots of the system
     # Note the difference between continuous and discrete-time systems...
-    minreal(SisoZpk(z,p,k), sqrt(eps()))
+    minreal(SisoZpk(z::Vector{TRnew},p::Vector{TRnew},k), sqrt(eps())) # type assert required or inference
 end
 
 

src/types/SisoTfTypes/conversion.jl

@@ -21,7 +21,7 @@ end
 
 
 function Base.convert(::Type{SisoRational{T}}, f::SisoRational) where {T<:Number}
-    return SisoRational{T}(convert(Poly{T}, f.num), convert(Poly{T}, f.den))
+    return SisoRational{T}(convert(Polynomial{T}, f.num), convert(Polynomial{T}, f.den))
 end
 
 function Base.convert(::Type{SisoRational{T1}}, f::SisoZpk{T2,TR}) where {T1<:Number, T2<:Number,TR<:Number}

src/types/TransferFunction.jl

@@ -113,7 +113,7 @@ function isapprox(G1::TransferFunction, G2::TransferFunction; kwargs...)
 end
 
 function isapprox(G1::Array{T}, G2::Array{S}; kwargs...) where {T<:SisoTf, S<:SisoTf}
-    all(i -> isapprox(G1[i], G2[i]; kwargs...), eachindex(G1))
+    all(i -> isapprox(promote(G1[i], G2[i])...; kwargs...), eachindex(G1))
 end
 
 ## ADDITION ##

src/types/conversion.jl

@@ -1,7 +1,7 @@
 # Base.convert(::Type{<:SisoTf}, b::Real) = Base.convert(SisoRational, b)
 # Base.convert{T<:Real}(::Type{<:SisoZpk}, b::T) = SisoZpk(T[], T[], b)
 # Base.convert{T<:Real}(::Type{<:SisoRational}, b::T) = SisoRational([b], [one(T)])
-# Base.convert{T1}(::Type{SisoRational{Vector{T1}}}, t::SisoRational) =  SisoRational(Poly(T1.(t.num.a)),Poly(T1.(t.den.a)))
+# Base.convert{T1}(::Type{SisoRational{Vector{T1}}}, t::SisoRational) =  SisoRational(Polynomial(T1.(t.num.coeffs$1),Polynomial(T1.(t.den.coeffs$1))
 # Base.convert(::Type{<:StateSpace}, t::Real) = ss(t)
 #
 
@@ -125,7 +125,7 @@ function siso_tf_to_ss(T::Type, f::SisoRational)
     # Get numerator coefficient of the same order as the denominator
     bN = length(num) == N+1 ? num[1] : 0
 
-    if N==0 #|| num == zero(Poly{T})
+    if N==0 #|| num == zero(Polynomial{T})
         A = fill(zero(T), 0, 0)
         B = fill(zero(T), 0, 1)
         C = fill(zero(T), 1, 0)
@@ -299,7 +299,7 @@ end
 #     λ = eigvals(A);
 #     T = promote_type(primitivetype(A), Float64)
 #     I = one(T)
-#     p = reduce(*,Poly([I]), Poly[Poly([I, -λᵢ]) for λᵢ in λ]);
+#     p = reduce(*,Polynomial([I]), Polynomial[Polynomial([I, -λᵢ]) for λᵢ in λ]);
 #     if maximum(imag.(p[:])./(I+abs.(real.(p[:])))) < sqrt(eps(T))
 #         for i = 1:length(p)
 #             p[i] = real(p[i])

src/types/promotion.jl

@@ -1,4 +1,4 @@
-# Base.eltype{T}(::Type{Poly{T}}) = T
+# Base.eltype{T}(::Type{Polynomial{T}}) = T
 # Base.eltype{T}(::Type{TransferFunction{T}}) = T
 # Base.eltype{T}(::Type{SisoZpk{T}}) = T
 # Base.eltype{T}(::Type{SisoRational{T}}) = T
@@ -144,5 +144,5 @@ Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1, T2<:Number} =
 # Base.zero(::SisoTf) = zero(SisoRational)
 # Base.zero(::Type{<:SisoZpk}) = SisoZpk(Float64[],Float64[],0.0)
 # Base.zero(::SisoZpk) = Base.zero(SisoZpk)
-# Base.zero{T}(::Type{SisoRational{T}}) = SisoRational(zero(Poly{T}), one(Poly{T})) # FIXME: Is the in analogy with how zero for SisoRational is createdzer?
+# Base.zero{T}(::Type{SisoRational{T}}) = SisoRational(zero(Polynomial{T}), one(Polynomial{T})) # FIXME: Is the in analogy with how zero for SisoRational is createdzer?
 # Base.zero{T}(::SisoRational{T}) = Base.zero(SisoRational{T})

src/types/tf.jl

@@ -71,7 +71,7 @@ function tf(var::AbstractString, Ts::Real)
 end
 
 ## Constructors for polynomial inputs
-function tf(num::AbstractArray{PT}, den::AbstractArray{PT},  Ts::Real=0.0) where {T<:Number, PT <: Polynomials.Poly{T}}
+function tf(num::AbstractArray{PT}, den::AbstractArray{PT},  Ts::Real=0.0) where {T<:Number, PT <: Polynomials.Polynomial{T}}
     ny, nu = size(num, 1), size(num, 2)
     if (ny, nu) != (size(den, 1), size(den, 2))
         error("num and den dimensions must match")
@@ -86,6 +86,6 @@ function tf(num::AbstractArray{PT}, den::AbstractArray{PT},  Ts::Real=0.0) where
     return TransferFunction{SisoRational{T}}(matrix, Ts)
 end
 
-function tf(num::PT, den::PT,  Ts::Real=0.0) where {T<:Number, PT <: Polynomials.Poly{T}}
+function tf(num::PT, den::PT,  Ts::Real=0.0) where {T<:Number, PT <: Polynomials.Polynomial{T}}
     tf(fill(num,1,1), fill(den,1,1), Ts)
 end

src/utilities.jl

@@ -33,7 +33,7 @@ end
 """ f = printpolyfun(var)
 `fun` Prints polynomial in descending order, with variable `var`
 """
-printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printpoly(io, p, mimetype, descending_powers=true, var=var)
+printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printpoly(io, Polynomial(p.coeffs, var), mimetype, descending_powers=true)
 
 # NOTE: Tolerances for checking real-ness removed, shouldn't happen from LAPACK?
 # TODO: This doesn't play too well with dual numbers..
@@ -44,12 +44,12 @@ printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printp
 # returned by LAPACK routines for eigenvalues.
 function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
     T = real(cT)
-    poly_factors = Vector{Poly{T}}()
+    poly_factors = Vector{Polynomial{T}}()
     for k=1:length(roots)
         r = roots[k]
 
         if isreal(r)
-            push!(poly_factors,Poly{T}([-real(r),1]))
+            push!(poly_factors,Polynomial{T}([-real(r),1]))
         else
             if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
                 continue
@@ -59,7 +59,7 @@ function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
                 throw(AssertionError("Found pole without matching conjugate."))
             end
 
-            push!(poly_factors,Poly{T}([real(r)^2+imag(r)^2, -2*real(r), 1]))
+            push!(poly_factors,Polynomial{T}([real(r)^2+imag(r)^2, -2*real(r), 1]))
             # k += 1 # Skip one iteration in the loop
         end
     end
@@ -68,7 +68,7 @@ function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
 end
 # This function should hande both Complex as well as symbolic types
 function roots2poly_factors(roots::Vector{T}) where T <: Number
-    return [Poly{T}([-r, 1]) for r in roots]
+    return [Polynomial{T}([-r, 1]) for r in roots]
 end
 
 
@@ -99,7 +99,7 @@ function _fix_conjugate_pairs!(v::AbstractVector{<:Real})
 end
 
 # Should probably try to get rif of this function...
-poly2vec(p::Poly) = p.a[1:end]
+poly2vec(p::Polynomial) = p.coeffs[1:end]
 
 
 function unwrap!(M::Array, dim=1)

test/runtests.jl

@@ -9,7 +9,6 @@ include("framework.jl")
 eye_(n) = Matrix{Int64}(I, n, n)
 
 my_tests = [
-            "test_delayed_systems",
             "test_statespace",
             "test_transferfunction",
             "test_zpk",
@@ -27,6 +26,7 @@ my_tests = [
             "test_lqg",
             "test_synthesis",
             "test_partitioned_statespace",
+            "test_delayed_systems",
             "test_demo_systems",
             ]
 

test/test_delayed_systems.jl

@@ -177,7 +177,7 @@ function y_expected(t, K)
       end
 end
 
-@test ystep ≈ y_expected.(t, K) atol = 1e-13
+@test ystep ≈ y_expected.(t, K) atol = 1e-10
 
 function dy_expected(t, K)
       if t < 1

test/test_synthesis.jl

@@ -12,7 +12,7 @@ S = [1]
 T = [1]
 @testset "minreal + feedback" begin
 @test isapprox(minreal(feedback(P,C),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5)
-@test isapprox(numpoly(minreal(feedback(L),1e-5))[1].a, numpoly(tf(1,[1,1]))[1].a)# This test is ugly, but numerical stability is poor for minreal
+@test isapprox(numpoly(minreal(feedback(L),1e-5))[1].coeffs, numpoly(tf(1,[1,1]))[1].coeffs)# This test is ugly, but numerical stability is poor for minreal
 @test feedback2dof(B,A,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
 @test feedback2dof(P,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
 @test isapprox(pole(minreal(tf(feedback(Lsys)),1e-5)) , pole(minreal(feedback(L),1e-5)), atol=1e-5)
@@ -22,7 +22,7 @@ Cint = tf([1,1],[1,0])
 Lint = P*C
 
 @test isapprox(minreal(feedback(Pint,Cint),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5) # TODO consider keeping minreal of Int system Int
-@test isapprox(numpoly(minreal(feedback(Lint),1e-5))[1].a, numpoly(tf(1,[1,1]))[1].a)# This test is ugly, but numerical stability is poor for minreal
+@test isapprox(numpoly(minreal(feedback(Lint),1e-5))[1].coeffs, numpoly(tf(1,[1,1]))[1].coeffs)# This test is ugly, but numerical stability is poor for minreal
 @test isapprox(pole(minreal(tf(feedback(Lsys)),1e-5)) , pole(minreal(feedback(L),1e-5)), atol=1e-5)
 
 

test/test_transferfunction.jl

@@ -142,9 +142,9 @@ D_diffTs = tf([1], [2], 0.1)
 @test_throws ErrorException [z 0]                     # Sampling time mismatch (inferec could be implemented)
 
 # Test polynomial inputs
-Poly = ControlSystems.Polynomials.Poly
+Polynomial = ControlSystems.Polynomials.Polynomial
 v1, v2, v3, v4 = [1,2,3], [4,5,6], [7,8], [9,10,11]
-p1, p2, p3, p4 = Poly.(reverse.((v1,v2,v3,v4)))
+p1, p2, p3, p4 = Polynomial.(reverse.((v1,v2,v3,v4)))
 
 S1v = tf(v1, v2)
 S2v = tf(v3, v4)