Type piracy in `hcat, vcat`

[baggepinnen]

Hello and thanks for this package!

I've found that this package commits type piracy on the functions hcat and vcat. The offending definitions are located here https://github.com/andreasvarga/DescriptorSystems.jl/blob/2745eda6a24a2334d82a81fb9c4247438cf4473d/src/connections.jl#L205

It's generally not a good idea to define a method for a function you do not own, with a signature containing only types you do not own. In this case, the functions are owned by Base, and this package defines methods with signatures containing only Base types, e.g., UniformScaling. This break functionality in completely unrelated packages as soon as DescriptorSystems is loaded.

More information on this subject is available here https://docs.julialang.org/en/v1/manual/style-guide/#Avoid-type-piracy

In my instance, loading DescriptorSystems breaks Symbolics.jl due to the call to vcat here

function TermInterface.arguments(a::Mul)
    a.sorted_args_cache[] !== nothing && return a.sorted_args_cache[]
    args = sort!([unstable_pow(k, v) for (k,v) in a.dict], lt=<ₑ)
    a.sorted_args_cache[] = isone(a.coeff) ? args : vcat(a.coeff, args)
end

dispatches to a broken method in DescriptorSystems instead of the correct method in Base.

[andreasvarga]

First of all thanks and sorry for this. I must confess that I had huge difficulties to make the concatenation functions running correctly, since myself I encountered similar issues (probably related to the problematic raised by you), where hcat, vcat and hvcat interferred with functions with the same names in packages like SparseArrays, LinearOperators, to mention only two of them. It took me a lot of time to arrive to a working version, where all input signatures have been correctly handled in my test (possibly other tests could still fail). My tuning approach was more a trial-and-error procedure rather than a deep understanding of the underlying issue. This is why, I am now even in a worse situation, since simply I forgot a lot of issues I encountered that time.

The aboove mentioned function tries to extend a similar function in LinearAlgebra/UniformScaling.jl . I wonder if the situation there is fundamentally different, because the hcat and vcat functions are redefined in similar way. Looking again to them, I observe that there are now new implementations, which somehow try to provide a hierarchy of argument calls (I am just guessing). I could try to mimic this approach, but I am not sure if this will lead to the expected cleaning.

Since you certainly have more experience handling these issue, I wonder if you have a proposal how to approach this issue. It would be also useful if you could provide the commands which lead to the failure, to be able to make tests with them.

[baggepinnen]

I think key is to never define a method with a signature that does not contain at least one of your own types. In the examples here

$f(A::Union{DescriptorStateSpace,AbstractNumOrArray,UniformScaling}...) 

this matches vcat with only UniformScaling which does not belong to this package. One solution is to instead define the method

$f(A::DescriptorStateSpace, B::Union{DescriptorStateSpace,AbstractNumOrArray,UniformScaling}...)

this is of course limited since it only matches if there is a DescriptorStateSpace as the first argument. One could extend this indefinitely by adding more methods, but it's a rather unsatisfying approach.

A perhaps better approach is to define how DescriptorStateSpace and AbstractNumOrArray are supposed to promote. I.e., promote(A::DescriptorStateSpace, B::AbstractNumOrArray) should return two DescriptorStateSpace objects. You accomplish this by defining the appropriate promote_rule and convert methods according to https://docs.julialang.org/en/v1/manual/conversion-and-promotion/ This way, you shouldn't have to define any other methods than vcat with only DescriptorStateSpace, since vcat with mixed types will call promote on the types, after which they are all DescriptorStateSpace.

[andreasvarga]

Could you send me an example with the mentioned failure?

[baggepinnen]

Here's a reproducing example

using ModelingToolkit
function testsys(; k, Ti=false, Td=false, wp=1, wd=1,
    Ni    = Ti == 0 ? Inf : √(Td / Ti),
    Nd    = 12,
    y_max = Inf,
    y_min = y_max > 0 ? -y_max : -Inf,
    name
)
    @parameters t
    @variables x(t)=0 u_r(t)=0 [input=true] u_y(t)=0 [input=true] y(t) [output=true] e(t)=0 ep(t)=0 ed(t)=0 ea(t)=0

    @parameters k=k Td=Td wp=wp wd=wd Ni=Ni Nd=Nd
    eqs = [
        e ~ u_r - u_y 
        ep ~ wp*u_r - u_y  
        ed ~ wd*u_r - u_y  
    ]
    ODESystem(eqs, t, name=name)
end

testsys(k=1, Ti=1, Td = 1, name=:test) # works
using DescriptorSystems
testsys(k=1, Ti=1, Td = 1, name=:test) # broken

with the output

julia> testsys(k=1, Ti=1, Td = 1, name=:test)
ERROR: MethodError: no constructors have been defined for Any
Stacktrace:
  [1] to_matrix(T::Type{Any}, A::Int64, wide::Bool) (repeats 2 times)
    @ DescriptorSystems ~/.julia/packages/DescriptorSystems/UKVXU/src/dsutils.jl:5
  [2] promote_to_system_(n::Int64, #unused#::Type{Any}, A::Int64, Ts::Int64)
    @ DescriptorSystems ~/.julia/packages/DescriptorSystems/UKVXU/src/connections.jl:310
  [3] promote_to_systems(Ts::Int64, n::Vector{Int64}, k::Int64, #unused#::Type{Any}, A::Int64, B::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}})
    @ DescriptorSystems ~/.julia/packages/DescriptorSystems/UKVXU/src/connections.jl:314
  [4] vcat(::Int64, ::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}})
    @ DescriptorSystems ~/.julia/packages/DescriptorSystems/UKVXU/src/connections.jl:227
  [5] arguments(a::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/5dtUX/src/types.jl:808
  [6] vars!(vars::Set{Any}, O::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}; op::Type)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:264
  [7] vars!(vars::Set{Any}, O::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}; op::Type)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:265
  [8] #38
    @ ~/.julia/dev/ModelingToolkit/src/utils.jl:247 [inlined]
  [9] BottomRF
    @ ./reduce.jl:81 [inlined]
 [10] _foldl_impl(op::Base.BottomRF{ModelingToolkit.var"#38#39"{DataType}}, init::Set{Any}, itr::Vector{SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}})
    @ Base ./reduce.jl:58
 [11] foldl_impl
    @ ./reduce.jl:48 [inlined]
 [12] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
 [13] #mapfoldl#214
    @ ./reduce.jl:160 [inlined]
 [14] #foldl#215
    @ ./reduce.jl:178 [inlined]
 [15] #vars#37
    @ ~/.julia/dev/ModelingToolkit/src/utils.jl:247 [inlined]
 [16] #vars#36
    @ ~/.julia/dev/ModelingToolkit/src/utils.jl:246 [inlined]
 [17] vars
    @ ~/.julia/dev/ModelingToolkit/src/utils.jl:246 [inlined]
 [18] collect_vars!(states::OrderedCollections.OrderedSet{Any}, parameters::OrderedCollections.OrderedSet{Any}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:307
 [19] ODESystem(eqs::Vector{Equation}, iv::Num; kwargs::Base.Iterators.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:name,), Tuple{Symbol}}})
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/diffeqs/odesystem.jl:172
 [20] testsys(; k::Int64, Ti::Int64, Td::Int64, wp::Int64, wd::Int64, Ni::Float64, Nd::Int64, y_max::Float64, y_min::Float64, name::Symbol)
    @ Main ./REPL[6]:17
 [21] top-level scope
    @ REPL[9]:1

[andreasvarga]

Thanks. I think the following reproduces the above error without loading additional packages:

julia> vcat(1,[1;'s'])
3-element Vector{Any}:
 1
 1
  's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

julia> using DescriptorSystems
[ Info: Precompiling DescriptorSystems [a81e2ce2-54d1-11eb-2c75-db236b00f339]

julia> vcat(1,[1;'s'])
ERROR: MethodError: no constructors have been defined for Any
Stacktrace:
 [1] to_matrix(T::Type{Any}, A::Int64, wide::Bool) (repeats 2 times)
   @ DescriptorSystems ~\Documents\software\Julia\DescriptorSystems.jl\src\dsutils.jl:5
 [2] promote_to_system_(n::Int64, #unused#::Type{Any}, A::Int64, Ts::Int64)
   @ DescriptorSystems ~\Documents\software\Julia\DescriptorSystems.jl\src\connections.jl:302
 [3] promote_to_systems(Ts::Int64, n::Vector{Int64}, k::Int64, #unused#::Type{Any}, A::Int64, B::Vector{Any})
   @ DescriptorSystems ~\Documents\software\Julia\DescriptorSystems.jl\src\connections.jl:306
 [4] vcat(::Int64, ::Vector{Any})
   @ DescriptorSystems ~\Documents\software\Julia\DescriptorSystems.jl\src\connections.jl:219
 [5] top-level scope
   @ REPL[2]:1

I will try to fix this case by restricting vcat/hcat to numerical data. This fix should also cover your case.

Unfortunately, I am afraid the problematic of type piracy raised by you still remains. I will try to invest some time to get rid of it.

I will submit a new release for DescriptorSystems next week. I would appreciate having a feedback from you on this issue.

[andreasvarga]

I made a new version, which you can test. I hope it will work for you with other packages.

My way to alleviate type piracy is to catch the relevant cases (e.g., by testing there is no DescriptorSystemState object among the input parameters) and then call LinearAlgebra tools to perform hcat/vcat/hvcat. You would probably not recommend this approach, but it does the work just now. I am still reflecting how to implement the promotion-based approach you suggested (actually this is already implemented in DescriptorSystems for handling hcat/vcat/hvcat-functionality for rational transfer functions). I have probably to change the definition of DescriptorSystemState object to include the sampling time in the type definition.

During my tests, I had to face the case [1 2; I] which normally should work with LinearAlgebra alone.

Surprisingly, the following crashes

julia> using LinearAlgebra

julia> [1 2; I]
ERROR: ArgumentError: number of columns of each array must match (got (2, 1))
Stacktrace:
 [1] _typed_vcat(#unused#::Type{Any}, A::Tuple{Matrix{Any}, Matrix{Any}})
   @ Base .\abstractarray.jl:1553
 [2] typed_vcat(::Type{Any}, ::Matrix{Any}, ::Matrix{Any})
   @ Base .\abstractarray.jl:1567
 [3] typed_hvcat(::Type{Any}, ::Tuple{Int64, Int64}, ::Int64, ::Vararg{Any, N} where N)
   @ Base .\abstractarray.jl:1957
 [4] hvcat(::Tuple{Int64, Int64}, ::Int64, ::Vararg{Any, N} where N)
   @ Base .\abstractarray.jl:1932
 [5] top-level scope
   @ REPL[9]:1

However, with the DescriptorSystems's type piracy, it works correctly: julia> using DescriptorSystems [ Info: Precompiling DescriptorSystems [a81e2ce2-54d1-11eb-2c75-db236b00f339]

julia> [1 2; I]
3×2 Matrix{Int64}:
 1  2
 1  0
 0  1

[baggepinnen]

I have probably to change the definition of DescriptorSystemState object to include the sampling time in the type definition.

In ControlSystems.jl, we have gone the route of including the time-domain in the type information, i.e., Continuous/Discrete is part of the type, but the numeric value of the same time is not. It's important to not include the sample time value itself, since then you would incur compilation of all methods for each unique sample time. The only time our approach is inadequate is in situations when you'd like to promote the I in [I sysd], since you do not know which sample time to give the static system I due to the sample time not being part of the type. One solution to this problem would be to allow the sample time -1 to denote "inherited" sample time, i.e., that the systems that have sample time -1 should adopt the sample times of the surrounding systems. I ave yet to implement the "inherited" sample time in ControlSystems.jl, but for all other purposes, our model has worked out well.

[andreasvarga]

I am reflecting to use instead the unspecified sample time -1, the value 1. This covers both the case of a sample time Ts explicitly specified as 1, but also the case of unspecified sampling time, which would correspond to a discrete-time system of the formx(t+1) = Ax(t)+Bu(t). Note that in frequency response computations (also in Matlab), this is implicitly used since abs(Ts) is always used instead just Ts. Thus the sample time -1 could be used in the sense you proposed only for inheritance purposes.

[dkarrasch]

To fix these kind of issues (I've struggled with that a lot in LinearMaps.jl), you should look up the corresponding method in LinearAlgebra and then add your own type(s) to the list in the union. Like this

# LinearAlgebra/uniformscaling.jl
for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{AbstractVecOrMat,UniformScaling}...) = ...
        @inline $f(A::Union{AbstractVecOrMat,UniformScaling,Number}...) = ...

Then you can define

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorSystem,AbstractVecOrMat,UniformScaling,Number}...) = # your_own_code

Since the union is strictly larger than the union in LinearAlgebra, your method is going to be more generic/less specific. That means, that your method gets called only in the case when there is at least one ::DescriptorSystem object in the list of arguments. From a quick glance, you seem to do it right, but it may help to take the copy-paste-(add-to-the-union) approach (instead of using your own const abbreviation). You could then check which method gets called via @which hcat(A, B, C) for different types, including and excluding objects of your own type. That's at least how I finally solved the issue after a painful tour in LinearMaps.jl.

[andreasvarga]

Thanks. Yes, this is the way I will try to address the issue.

[andreasvarga]

In Julia 1.7 the code below is not yet included. Should I wait for Julia 1.8 with the fix?

Daniel Karrasch @.***> schrieb am Fr., 4. März 2022, 10:23:

To fix these kind of issues (I've struggled with that a lot in LinearMaps.jl), you should look up the corresponding method in LinearAlgebra and then add your own type(s) to the list in the union. Like this

LinearAlgebra/uniformscaling.jlfor (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))

@eval begin
    @inline $f(A::Union{AbstractVecOrMat,UniformScaling}...) = ...
    @inline $f(A::Union{AbstractVecOrMat,UniformScaling,Number}...) = ...

Then you can define

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols")) @eval begin @inline $f(A::Union{DescriptorSystem,AbstractVecOrMat,UniformScaling,Number}...) = # your_own_code

Since the union is strictly larger than the union in LinearAlgebra, your method is going to be more generic/less specific. That means, that your method gets called only in the case when there is at least one ::DescriptorSystem object in the list of arguments. From a quick glance, you seem to do it right, but it may help to take the copy-paste-(add-to-the-union) approach (instead of using your own const abbreviation). You could then check which method gets called via @which hcat(A, B, C) for different types, including and excluding objects of your own type. That's at least how I finally solved the issue after a painful tour in LinearMaps.jl.

— Reply to this email directly, view it on GitHub https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1058987339, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEE5OC6TVLIRRQHMBV3U6HJB5ANCNFSM5HLI6N5Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[dkarrasch]

The question is if you need to support concatenation with numbers. If not, you could simply leave out Number above and proceed as described. I would actually do that to see if that fixes the issues. Supporting Numbers makes sense perhaps only when the stdlib supports it, so starting from v1.8.

[andreasvarga]

This is what I implemented following your suggestions to avoid type piracy (probably in a wrong way, see the results):

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace, Number}...) = $_f(A...)
        @inline $f(A::Union{DescriptorStateSpace, UniformScaling}...) = $_f(A...)
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat{<:Number}, UniformScaling}...) = $_f(A...)
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat{<:Number},UniformScaling, Number}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat{<:Number},UniformScaling, Number}...)
            n = -1
            for a in A
                if !isa(a, UniformScaling)
                    require_one_based_indexing(a)
                    na = size(a,$dim)
                    n >= 0 && n != na &&
                        throw(DimensionMismatch(string("number of ", $name,
                            " of each array must match (got ", n, " and ", na, ")")))
                    n = na
                end
            end
            n == -1 && throw(ArgumentError($("$f of only UniformScaling objects cannot determine the matrix size")))
            #if isadss(A...) == 0
            if isnothing(findfirst(a -> isa(a, DescriptorStateSpace), A)) 
                println("type piracy")
                # alleviate type piracy problematic
                n == 1 && (A = promote_to_arrays(A...))  # convert all scalars to vectors
                return cat(LinearAlgebra.promote_to_arrays(fill(n, length(A)), 1, Matrix, A...)..., dims=Val(3-$dim))
            else       
                Ts = promote_system_SamplingTime(A...)
                return $f(promote_to_systems(Ts, fill(n,length(A)), 1, promote_type(eltype.(A)...), A...)...)
            end
        end
    end
end

I am sad to report the following disillusioning results:

julia> hcat(rand(1,2), rand(1,2))
1×4 Matrix{Float64}:
 0.721178  0.297783  0.853222  0.0502695

julia> hcat(rand(1,2), rand(1,2),1)
type piracy
1×5 Matrix{Float64}:
 0.210189  0.781186  0.0382634  0.385548  1.0

julia> hcat(rand(1,2), rand(1,2),I)
type piracy
1×5 Matrix{Float64}:
 0.86124  0.961221  0.326077  0.783298  1.0

julia> hcat(rand(1,2), 1, I)
type piracy
1×4 Matrix{Float64}:
 0.641297  0.128555  1.0  1.0

julia> hcat(1, 1, I)
type piracy
1×3 Matrix{Int64}:
 1  1  1

I hope you have a salvatory idea, what I did wrong.

[dkarrasch]

You need to have the exact same signature as the method in LinearAlgebra, and then add DescriptorSystem. By restricting the arrays to Number eltypes, you change the specificity and make your method more specific, and hence it takes over precedence.

[andreasvarga]

Since I am detecting and appropriately handling type piracy cases, is this still type piracy?

Daniel Karrasch @.***> schrieb am Fr., 4. März 2022, 22:41:

You need to have the exact same signature as the method in LinearAlgebra, and then add DescriptorSystem. By restricting the arrays to Number eltypes, you change the specificity and make your method more specific, and hence it takes over precedence.

— Reply to this email directly, view it on GitHub https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1059551906, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEFRSC6RB5O6GPZKWHDU6J7QDANCNFSM5HLI6N5Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[andreasvarga]

I changed to

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{AbstractVecOrMat, UniformScaling}...) = $_f(A...)
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...)
            n = -1
            for a in A
                if !isa(a, UniformScaling)
                    require_one_based_indexing(a)
                    na = size(a,$dim)
                    n >= 0 && n != na &&
                        throw(DimensionMismatch(string("number of ", $name,
                            " of each array must match (got ", n, " and ", na, ")")))
                    n = na
                end
            end
            n == -1 && throw(ArgumentError($("$f of only UniformScaling objects cannot determine the matrix size")))
            #if isadss(A...) == 0
            if isnothing(findfirst(a -> isa(a, DescriptorStateSpace), A)) 
                println("type piracy")
                # alleviate type piracy problematic
                n == 1 && (A = promote_to_arrays(A...))  # convert all scalars to vectors
                return cat(LinearAlgebra.promote_to_arrays(fill(n, length(A)), 1, Matrix, A...)..., dims=Val(3-$dim))
            else       
                Ts = promote_system_SamplingTime(A...)
                return $f(promote_to_systems(Ts, fill(n,length(A)), 1, promote_type(eltype.(A)...), A...)...)
            end
        end
    end
end

with practically the same results as before:

julia> hcat(rand(1,2), rand(1,2),I)
type piracy
1×5 Matrix{Float64}:
 0.946567  0.484623  0.870328  0.350247  1.0

julia> hcat(I,rand(1,2),[1])
type piracy
1×4 Matrix{Float64}:
 1.0  0.555072  0.427686  1.0

julia> hcat(I,rand(1,2),1)
1×4 Matrix{Any}:
 UniformScaling{Bool}(true)  0.590048  0.230466  1

julia> @which hcat(I,rand(1,2),1)
hcat(X...) in Base at abstractarray.jl:1823

The results in the first two cases should not occur, while the third case is not handled correctly in Julia 1.7!

Additionally I observed the following type piracy (?) (unrelated to DescriptorSystems):

julia> @which hcat(rand(1,2), rand(1,2),rand(1,2))
hcat(A::Union{Vector{T}, Matrix{T}, Adjoint{T, Vector{T}}, Transpose{T, Vector{T}}, LinearAlgebra.AbstractTriangular{T, A} where A<:(Matrix), Hermitian{T, A} where A<:(Matrix), Symmetric{T, A} where A<:(Matrix)}...) where T in SparseArrays at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\SparseArrays\src\sparsevector.jl:1131

[dkarrasch]

Sure, because you overwrite this method:

@inline $f(A::Union{AbstractVecOrMat, UniformScaling}...) = $_f(A...)

Again, you need to take one concrete method from LinearAlgebra and add DescriptorSystem to the union in the signature of that very method, to make sure that your method is less specific than the one you picked. The type piracy is that you redefine methods that this package doesn't own with signature consisting exclusively of types that this package doesn't own.

while the third case is not handled correctly in Julia 1.7!

Yes, this is fixed in v1.8. That was a longstanding issue.

Additionally I observed the following type piracy (?) (unrelated to DescriptorSystems):

And this type piracy has been resolved by decoupling SparseArrays from LinearAlgebra in v1.8.

[andreasvarga]

Again, you need to take one concrete method from LinearAlgebra and add DescriptorSystem to the union in the signature of that very method, to make sure that your method is less specific than the one you picked.

Would it be possible to show me how this can be done? I had the impression that I did exacly what you suggested in:

To fix these kind of issues (I've struggled with that a lot in LinearMaps.jl), you should look up the corresponding method in LinearAlgebra and then add your own type(s) to the list in the union. Like this

# LinearAlgebra/uniformscaling.jl
for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{AbstractVecOrMat,UniformScaling}...) = ...
        @inline $f(A::Union{AbstractVecOrMat,UniformScaling,Number}...) = ...

[dkarrasch]

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace,AbstractVecOrMat,UniformScaling}...) = ...
        @inline $f(A::Union{DescriptorStateSpace,AbstractVecOrMat,UniformScaling,Number}...) = ...

perhaps without the second method (or the second method only for VERSION>=v"1.8"; BTW: it is exactly this second method that fixes your "third case" from above).

[andreasvarga]

I tried only with

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace,AbstractVecOrMat,UniformScaling}...) = ...

(see the code above).

What I have to put instead of ... ? I wonder if I really need the additional function _f or I could try to manage without this function?

[dkarrasch]

Yes, you can put your code from _*cat directly there. The reason these two functions (which both referred to _*cat are there is exactly multiple dispatch and the fact that several packages extended $f(A::Union{AbstractVecOrMat,UniformScaling}...) and needed a function in LinearAlgebra to avoid ambiguity.

[andreasvarga]

Yes, you can put your code from _*cat directly there.

If you look to my code, I did exactly that.

[dkarrasch]

In https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1059547899, you added the {<:Number} restriction to AbstractVecOrMat, in https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1059729954 you redefined @inline $f(A::Union{AbstractVecOrMat, UniformScaling}...) = $_f(A...). That's why I suggested the copy-paste-add approach: not changing anything from the signature, and only add your own type to the union (of that method, and not define an additional method with your type added). You should only define one method as in https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1059741715 (and potentially another one with Numbers added), and you can put your actual code instead of ....

[andreasvarga]

I implemented the following interface, following your suggestions:

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{AbstractVecOrMat, UniformScaling}...) = $_f(A...)
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...)

Before giving up, I would like to recap what this setup should provide:

1. The command hcat(rand(1,2),I) should activate the LinearAlgebra functionhcat. Unfortunately, this doesn't happen:

   julia> @which hcat(rand(1,2),I)
   hcat(A::Union{UniformScaling, AbstractVecOrMat{T} where T}...) in DescriptorSystems at C:\Users\Andreas\Documents\software\Julia\DescriptorSystems.jl\src\connections.jl:234

2. The command hcat(rss(0,1,1), [1],1) should not work (the Number option not covered). This happens indeed

   
   julia> hcat(rss(0,1,1), [1],1)
   1×3 Matrix{DescriptorStateSpace{Float64}}:
   DescriptorStateSpace{Float64}

Feedthrough matrix D: 1×1 Matrix{Float64}: 0.809985

Static gain. … DescriptorStateSpace{Float64}

Feedthrough matrix D: 1×1 Matrix{Float64}: 1.0

Static gain.

which is a 1x3 `Matrix` of Descriptor Systems (instead of a single Descriptor System).

3. The command  
`hcat(rss(0,1,1), [1],I)`
should work. Indeed

julia> hcat(rss(0,1,1), [1],I) DescriptorStateSpace{Float64}

Feedthrough matrix D: 1×3 Matrix{Float64}: 0.802032 1.0 1.0

Static gain.

If you have other suggestions, I will be happy to check them.

[dkarrasch]

Don't implement

@inline $f(A::Union{AbstractVecOrMat, UniformScaling}...) = $_f(A...)

this overrides LinearAlgebra behavior, just implement

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling}...)

or put the $_f code in the $f code block right away.

[andreasvarga]

With the following interface

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...)

I obtained the following results:

julia> using LinearAlgebra

julia> vcat(1,[1;'s'])
3-element Vector{Any}:
 1
 1
  's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

julia> using DescriptorSystems
[ Info: Precompiling DescriptorSystems [a81e2ce2-54d1-11eb-2c75-db236b00f339]

julia> vcat(1,[1;'s'])
*** type piracy ***
3-element Vector{Any}:
 1
 1
  's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

This was the motivating case for this issue.

The good news:

julia> hcat(rand(1,2), rand(1,2))
1×4 Matrix{Float64}:
 0.833964  0.795188  0.0312341  0.548431

julia> hcat(rand(1,2), rand(1,2),I)
1×5 Matrix{Float64}:
 0.549049  0.504214  0.836734  0.831469  1.0

The not so good news (possibly will be handled in Julia 1.8 correctly):

julia> hcat(rand(1,2), rand(1,2),1)
*** type piracy ***
1×5 Matrix{Float64}:
 0.121579  0.758398  0.464527  0.981774  1.0

julia> hcat(rand(1,2), 1, I)
*** type piracy ***
1×4 Matrix{Float64}:
 0.334563  0.372151  1.0  1.0

julia> hcat(1, 1, I)
*** type piracy ***
1×3 Matrix{Int64}:
 1  1  1

[dkarrasch]

With the following interface

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))
    @eval begin
        @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...) = $_f(A...)
        function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...)

I obtained the following results:

julia> using LinearAlgebra

julia> vcat(1,[1;'s'])
3-element Vector{Any}:
 1
 1
  's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

julia> using DescriptorSystems
[ Info: Precompiling DescriptorSystems [a81e2ce2-54d1-11eb-2c75-db236b00f339]

julia> vcat(1,[1;'s'])
*** type piracy ***
3-element Vector{Any}:
 1
 1
  's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

This was the motivating case for this issue.

The not so good news (possibly will be handled in Julia 1.8 correctly):

julia> hcat(rand(1,2), rand(1,2),1)
*** type piracy ***
1×5 Matrix{Float64}:
 0.121579  0.758398  0.464527  0.981774  1.0

julia> hcat(rand(1,2), 1, I)
*** type piracy ***
1×4 Matrix{Float64}:
 0.334563  0.372151  1.0  1.0

julia> hcat(1, 1, I)
*** type piracy ***
1×3 Matrix{Int64}:
 1  1  1

These classes of piracy issues all arise because you included Number in the signature on v1.7 (against my recommendation 😉 ). This is not wise because that doesn't work properly on v1.7 (with UniformScalings) even in the absence of DescriptorStateSpace. If you remove that and avoid such test cases, then you should be fine. If you want the Number stuff to work, you (and everybody else) will need to upgrade to v1.8+ anyway.

[andreasvarga]

In a previous message I emphasized that I am detecting type piracy and perform the cat operations with appropriate linear algebra tools. Therefore, I think this is not a type piracy at all.

The functionality with Number is essential for the functional consitency of model building, so I need to keep it. I hope to have a clean version with 1.8.

An issue for which I don't have an answer is how to restrict all cat operations to numerical matrices and avoid generated constructs of type Any.

Daniel Karrasch @.***> schrieb am Sa., 5. März 2022, 19:38:

With the following interface

for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"))

@eval begin

    @inline $f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...) = $_f(A...)

    function $_f(A::Union{DescriptorStateSpace, AbstractVecOrMat,UniformScaling,Number}...)

I obtained the following results:

julia> using LinearAlgebra

julia> vcat(1,[1;'s'])

3-element Vector{Any}:

1

1

's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

julia> using DescriptorSystems

[ Info: Precompiling DescriptorSystems [a81e2ce2-54d1-11eb-2c75-db236b00f339]

julia> vcat(1,[1;'s'])

type piracy

3-element Vector{Any}:

1

1

's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

This was the motivating case for this issue.

The not so good news (possibly will be handled in Julia 1.8 correctly):

julia> hcat(rand(1,2), rand(1,2),1)

type piracy

1×5 Matrix{Float64}:

0.121579 0.758398 0.464527 0.981774 1.0

julia> hcat(rand(1,2), 1, I)

type piracy

1×4 Matrix{Float64}:

0.334563 0.372151 1.0 1.0

julia> hcat(1, 1, I)

type piracy

1×3 Matrix{Int64}:

1 1 1

These classes of piracy issues all arise because you included Number in the signature on v1.7 (against my recommendation 😉 ). This is not wise because that doesn't work properly on v1.7 (with UniformScalings) even in the absence of DescriptorStateSpace. If you remove that and avoid such test cases, then you should be fine. If you want the Number stuff to work, you (and everybody else) will need to upgrade to v1.8+ anyway.

— Reply to this email directly, view it on GitHub https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1059812578, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHECY2K5LC3LIVJITGWDU6OSZRANCNFSM5HLI6N5Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[dkarrasch]

That's fine. If you need it to work with matrices, uniform scalings, numbers and DescriptorStateSpace, then you should require Julia v1.8, since even the matrix-uniform scaling-number combination doesn't work in Julia v1.7. In that case, there is no need to catch the type piracy case, which it remains, because multiple dispatch redirects function calls, even though you seem to handle that, but who gives a guarantee, and how do you keep up with potential development of LinearAlgebra? Avoiding type piracy means exactly this: don't override behavior, don't redirect function calls in multiple dispatch by just loading a package, or otherwise you get issues like the original one (Symbolics.jl breaks because the vcat from DescriptorSystems.jl, which has nothing to do with Symbolics.jl, gets called).

An issue for which I don't have an answer is how to restrict all cat operations to numerical matrices and avoid generated constructs of type Any.

I don't think you should be concerned about that. If the user provides strange data, they get strange results. And who knows, maybe the result isn't strange. That's the point of generic programming: the user may use your functions that were intended for one case, but by method overload your code extends well to types that you have never thought of. So, if I understand correctly, as long as users do provide numeric matrices, they do get matrices with concrete numeric eltype with your concatenation function.

[andreasvarga]

I would like to conclude this productive discussion by thanking you for the help and support to address the type piracy issue caused by the usage of DescriptorSystems. Just for the record: to alleviate/fix this problem, I produced new versions of MatrixEquations, MatrixPencils and DescriptorSystems, which are free of any internal type piracy (mainly originating in MatrixPencils running with DescriptorSystems) in Julia 1.6, 1.7 and 1.8. The only type piracy occurs in specially constructed tests and manifests only under Julia 1.6 and 1.7 with an appropriate warning message. By including the handling of Number type, a side effect of this is that constructs, such as [1 2;I], which usually fail using only the LinearAlgebra, can be handled if DescriptorSystems is running

julia> [1 2;I]
┌ Warning: Type piracy in hvcat: to be fixed in Julia 1.8 (make an issue otherwise)
└ @ DescriptorSystems C:\Users\Andreas\Documents\software\Julia\DescriptorSystems.jl\src\connections.jl:306
3×2 Matrix{Int64}:
 1  2
 1  0
 0  1

This feature is useful as long as the new implementations of hcat/vcat/hvcat in Julia 1.8 are not backported to Julia 1.6 and 1.7.

In Julia 1.8, no such warnings are issued, so everything seems to be fine.

[andreasvarga]

I removed the warnings which became boring!