From 5230af1297904a6a6b73bafeb5b72af4b4b99242 Mon Sep 17 00:00:00 2001 From: Quentin Quadrat Date: Sat, 7 Sep 2019 03:09:52 +0200 Subject: [PATCH] WIP --- Manifest.toml | 49 ++++++ Project.toml | 1 + src/MaxPlus.jl | 396 +++++++++++++++++++++++++++++++++++------------ test/runtests.jl | 203 ++++++++++++++++++------ 4 files changed, 509 insertions(+), 140 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 081a570..a1ee87e 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,6 +1,20 @@ [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "2.1.0" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + [[Distributed]] deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -9,6 +23,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" deps = ["LinearAlgebra", "Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[LibGit2]] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -23,17 +40,35 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[Pkg]] +deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + [[Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + [[Random]] deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +[[SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -41,9 +76,23 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc" deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[Suppressor]] +deps = ["Compat"] +git-tree-sha1 = "a39342763981e766a72938b59032dc02a2d74da5" +uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +version = "0.1.1" + [[Test]] deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[UUIDs]] +deps = ["Random"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/Project.toml b/Project.toml index a4ddd28..fa4b15d 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.0" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" [compat] julia = "1" diff --git a/src/MaxPlus.jl b/src/MaxPlus.jl index 627a78e..28bd3c1 100644 --- a/src/MaxPlus.jl +++ b/src/MaxPlus.jl @@ -12,61 +12,214 @@ export MP, ⨁, ⨂, mpzero, mpone, mp0, mp1, mptop, mpI, mpeye, mpzeros, mpones, - mpsparse, full, dense, plustimes, array, mparray, mptrace, merde + mpsparse, full, dense, array, + plustimes, minplus, + mptrace, + mp_set_show_mode, LaTeX # ============================================================================== """ MP{T} -Max-Plus immutable structure. -Promote a number of type T (such as ideally Float64) to a -max-plus number. +Max-Plus immutable structure. Promote a number of type T (ideally a Float64) to +a max-plus number. # Examples ```julia-repl -a=MP(3) -typeof(a) +julia> a=MP(3.0); typeof(a) +MP{Float64} + +julia> a=MP(3); typeof(a) +MP{Int64} ``` """ struct MP{T} <: Number λ::T end # ============================================================================== - # Copy constructor + MP(x::MP) = MP(x.λ) # ============================================================================== +Base.promote_rule(::Type{MP{T}}, ::Type{U}) where T where U = MP{T} +Base.convert(::MP{T}, x) where T = MP(T(x)) + +# ============================================================================== + +""" + MP(A::Array) + +Convert a dense array to a max-plus dense array. + +# Examples +```julia-repl +julia> MP([1.0 -Inf; 0.0 4]) +2×2 Array{MP{Float64},2}: + 1.0 -Inf + 0.0 4.0 +``` +""" +MP(A::Array) = map(MP, A) + +# ============================================================================== + +""" + MP(A::SparseArray) + +Convert a sparse array to a max-plus sparse array. +By default values are max-plus zeros values are preserved else if +parameter `preserve` is set then + +# Examples +```julia-repl +julia> MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0])) + [1, 1] = -Inf + [2, 2] = 2.0 + [3, 3] = 0.0 + +julia> MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]), preserve=false) + [2, 2] = 2.0 +``` +""" +function MP(S::SparseMatrixCSC{T,U}; preserve=true) where T where U + if preserve + convert(SparseMatrixCSC{MP{T},U}, S) + else + M = spzeros(MP{T}, size(S,1), size(S,2)) + for c = 1:size(S, 2), r = nzrange(S, c) + if (S[r,c] != zero(T)) && (MP(S[r,c]) != mpzero(T)) + M[r,c] = convert(MP{T}, S[r,c]) + end + end + M + end +end + +# ============================================================================== + +""" + MP(A::SparseVector) + +Convert a sparse vector to a max-plus sparse vector. +By default values are max-plus zeros values are preserved else if +parameter `preserve` is set then + +# Examples +```julia-repl +julia> MP(sparse([1.0, 0.0, 1.0])) +3-element SparseVector{MP{Float64},Int64} with 2 stored entries: + [1] = 1.0 + [3] = 1.0 +``` +""" +function MP(V::SparseVector{T,U}; preserve=true) where T where U + if preserve + convert(SparseVector{MP{T},U}, V) + else + S = spzeros(MP{T}, size(V,1), size(V,2)) + for c = V.n + if (S[c] != zero(T)) && (MP(S[c]) != mpzero(T)) + M[c] = convert(MP{T}, S[c]) + end + end + S + end +end + +# ============================================================================== + +global mp_show_mode = 2 + +""" + mp_set_show_mode(mode::Int) + +Change the behavior of show(io::IO,x::MP). See help of show(io::IO,x::MP) +for examples. """ - show(io::IO,a::MP) +mp_set_show_mode(mode::Int) = (global mp_show_mode = min(max(mode, 0), 2)) + +""" + show(io::IO,x::MP) Display a max-plus number. -Inf are displayed with 'ε' symbols. Note in Scicoslab -Inf are displayed with '.' symbols. # Examples ```julia-repl -julia> mpeye(Float64, 5,5) -5×5 Array{MP{Float64},2}: - 0.000 ε ε ε ε - ε 0.000 ε ε ε - ε ε 0.000 ε ε - ε ε ε 0.000 ε - ε ε ε ε 0.000 +julia> mp_set_show_mode(0); mpeye(Float64, 2,2) +2×2 Array{MP{Float64},2}: + 0.0 -Inf + -Inf 0.0 + +julia> mp_set_show_mode(1); mpeye(Float64, 2,2) +2×2 Array{MP{Float64},2}: + 0.0 ε + ε 0.0 + +julia> mp_set_show_mode(2); mpeye(Float64, 2,2) +2×2 Array{MP{Float64},2}: + e ε + ε e ``` """ -Base.show(io::IO,a::MP) = - (a.λ == -Inf) ? (@printf io "ε") : show(io, a.λ) +function Base.show(io::IO, x::MP) + if mp_show_mode == 0 + show(io, x.λ) + elseif mp_show_mode == 1 + (x.λ == -Inf) ? (@printf io "ε") : show(io, x.λ) + else + (x.λ == -Inf) ? (@printf io "ε") : ((x.λ == 0) ? (@printf io "e") : show(io, x.λ)) + end +end + +# ============================================================================== + +""" + LaTeX(io::IO, A::Array{MP}) -#Base.show(io::IO,a::MP) = -# (a.λ == -Inf) ? (@printf io "ε") : ((a.λ == 0) ? (@printf io "e") : show(io, a.λ)) +Convert a max-plus dense matrix to a LaTeX formula + +# Examples +```julia-repl +julia> LaTeX(stdout, MP([4 3; 7 -Inf])) +\\left[ +\\begin{array}{*{20}c} +4 & 7 \\\\ +3 & \\varepsilon \\\\ +\\end{array} +\\right] +``` +""" +function LaTeX(io::IO, A::Array{MP{T}}) where T + (@printf io "\\left[\n\\begin{array}{*{20}c}\n") + for j = 1:size(A,2) + for i = 1:size(A,1) + if A[i,j] == mpzero(T) + (@printf io "\\varepsilon") + elseif A[i,j] == mpone(T) + (@printf io "e") + elseif (A[i,j].λ == trunc(A[i,j].λ)) + (@printf io "%d" A[i,j].λ) + else + (@printf io "%.3f" A[i,j].λ) + end + if i < size(A, 1) + (@printf io " & ") + end + end + (@printf io " \\\\\n") + end + (@printf io "\\end{array}\n\\right]\n") +end # ============================================================================== """ - +(a::MP, b::MP) + +(x::MP, y::MP) -Max operator. Return the maximum of `a` and `b`. +Max operator. Return the maximum of `x` and `y`. See also ⨁ operator. # Examples @@ -75,16 +228,16 @@ julia> MP(1.0) + MP(3.0) MP{Float64}(3.0) ``` """ -Base.:+(a::MP, b::MP) = MP(max(a.λ, b.λ)) -Base.:+(a::MP, b::Real) = MP(max(a.λ, b)) -Base.:+(a::Real, b::MP) = MP(max(a, b.λ)) +Base.:(+)(x::MP, y::MP) = MP(max(x.λ, y.λ)) +Base.:(+)(x::MP, y::Real) = MP(max(x.λ, y)) +Base.:(+)(x::Real, y::MP) = MP(max(x, y.λ)) # ============================================================================== """ - ⨁(a::MP, b::MP) + ⨁(x::MP, y::MP) -Max operator. Return the maximum of `a` and `b`. +Max operator. Return the maximum of `x` and `y`. Unicode character "circled plus": U+2A01. # Examples @@ -93,17 +246,17 @@ julia> ⨁(1.0, 3.0) MP{Float64}(3.0) ``` """ -⨁(a::MP, b::MP) = MP(max(a.λ, b.λ)) -⨁(a::MP, b::Real) = MP(max(a.λ, b)) -⨁(a::Real, b::MP) = MP(max(a, b.λ)) -⨁(a::Real, b::Real) = MP(max(a, b)) +⨁(x::MP, y::MP) = MP(max(x.λ, y.λ)) +⨁(x::MP, y::Real) = MP(max(x.λ, y)) +⨁(x::Real, y::MP) = MP(max(x, y.λ)) +⨁(x::Real, y::Real) = MP(max(x, y)) # ============================================================================== """ - *(a::MP, b::MP) + *(x::MP, y::MP) -Addition operator. Return the summation of `a` and `b`. +Addition operator. Return the summation of `x` and `y`. # Examples ```julia-repl @@ -111,16 +264,22 @@ julia> MP(1.0) * MP(3.0) MP{Float64}(4.0) ``` """ -Base.:*(a::MP, b::MP) = MP(a.λ + b.λ) -Base.:*(a::MP, b::Real) = MP(a.λ + b) -Base.:*(a::Real, b::MP) = MP(a + b.λ) +Base.:(*)(x::MP, y::MP) = MP(x.λ + y.λ) +Base.:(*)(x::MP, y::Real) = MP(x.λ + y) +Base.:(*)(x::Real, y::MP) = MP(x + y.λ) + +# ============================================================================== + +@inline Base.literal_pow(::typeof(^), x::MP, ::Val{0}) = one(x) +@inline Base.literal_pow(::typeof(^), x::MP, ::Val{p}) where {p} = MP(x.λ * p) +@inline Base.:(^)(x::MP, y::Number) = MP(x.λ * y) # ============================================================================== """ - ⨂(a::MP, b::MP) + ⨂(x::MP, y::MP) -Return the summation of `a` and `b`. +Return the summation of `x` and `y`. Unicode character "circled times": U+2A00. # Examples @@ -129,17 +288,17 @@ julia> ⨂(1.0, 3.0) MP{Float64}(4.0) ``` """ -⨂(a::MP, b::MP) = MP(a.λ + b.λ) -⨂(a::MP, b::Real) = MP(a.λ + b) -⨂(a::Real, b::MP) = MP(a + b.λ) -⨂(a::Real, b::Real) = MP(a + b) +⨂(x::MP, y::MP) = MP(x.λ + y.λ) +⨂(x::MP, y::Real) = MP(x.λ + y) +⨂(x::Real, y::MP) = MP(x + y.λ) +⨂(x::Real, y::Real) = MP(x + y) # ============================================================================== """ - -(a::MP, b::MP) + -(x::MP, y::MP) -Minus operator. Return the difference between `a` and `b`. +Minus operator. Return the difference between `x` and `y`. # Examples ```julia-repl @@ -147,16 +306,16 @@ julia> MP(1.0) / MP(2.0) MP{Float64}(-1.0) ``` """ -Base.:-(a::MP, b::MP) = MP(a.λ - b.λ) -Base.:-(a::MP, b::Real) = MP(a.λ - b) -Base.:-(a::Real, b::MP) = MP(a - b.λ) +Base.:(-)(x::MP, y::MP) = MP(x.λ - y.λ) +Base.:(-)(x::MP, y::Real) = MP(x.λ - y) +Base.:(-)(x::Real, y::MP) = MP(x - y.λ) # ============================================================================== """ - /(a::MP, b::MP) + /(x::MP, y::MP) -Divisor operator. Return the difference between `a` and `b`. +Divisor operator. Return the difference between `x` and `y`. # Examples ```julia-repl @@ -164,16 +323,16 @@ julia> MP(1.0) / MP(2.0) MP{Float64}(-1.0) ``` """ -Base.:/(a::MP, b::MP) = MP(a.λ - b.λ) -Base.:/(a::MP, b::Real) = MP(a.λ - b) -Base.:/(a::Real, b::MP) = MP(a - b.λ) +Base.:(/)(x::MP, y::MP) = MP(x.λ - y.λ) +Base.:(/)(x::MP, y::Real) = MP(x.λ - b) +Base.:(/)(x::Real, y::MP) = MP(a - y.λ) # ============================================================================== """ - min(a::MP, b::MP) + min(x::MP, y::MP) -Return the minimun of `a` and `b`. +Return the minimun of `x` and `y`. # Examples ```julia-repl @@ -181,15 +340,18 @@ julia> min(MP(1), -3) MP{Int64}(-3) ``` """ -Base.:min(a::MP, b::MP) = MP(min(a.λ, b.λ)) -Base.:min(a::MP, b::Real) = MP(min(a.λ, b)) -Base.:min(a::Real, b::MP) = MP(min(a, b.λ)) +Base.:min(x::MP, y::MP) = MP(min(x.λ, y.λ)) +Base.:min(x::MP, y::Real) = MP(min(x.λ, y)) +Base.:min(x::Real, y::MP) = MP(min(x, y.λ)) # ============================================================================== -Base.isless(a::MP{T}, b::MP{T}) where T = a.λ < b.λ -Base.promote_rule(::Type{MP{T}}, ::Type{U}) where T where U = MP{T} -Base.convert(::MP{T}, x) where T = MP(T(x)) +Base.:(==)(x::MP, y::Real) = (x.λ == y) +Base.:(==)(x::Real, y::MP) = (x == y.λ) + +Base.isless(x::MP, y::MP) = x.λ < y.λ +Base.isless(x::MP, y::Real) = x.λ < y +Base.isless(x::Real, y::MP) = x < y.λ # ============================================================================== @@ -199,8 +361,8 @@ Base.convert(::MP{T}, x) where T = MP(T(x)) Create the max-plus constant zero equal to -Inf. This value is the neutral for the ⨁ operator. """ -Base.zero(::MP{T}) where T = MP(typemin(T)) Base.zero(::Type{MP{T}}) where T = MP(typemin(T)) +Base.zero(x::MP{T}) where T = zero(typeof(x)) """ mpzero(::Type{T}) @@ -216,8 +378,8 @@ mpzero(::Type{T}) where T = MP(typemin(T)) Create the max-plus constant one equal to 0. This value is the neutral for operators ⨁ and ⨂. """ -Base.one(::MP{T}) where T = MP(zero(T)) Base.one(::Type{MP{T}}) where T = MP(zero(T)) +Base.one(x::MP{T}) where T = one(typeof(x)) """ one(::Type{T}) @@ -246,7 +408,7 @@ julia> mp0 + 5 5.0 ``` """ -mp0 = mpzero(Float64) +const global mp0 = mpzero(Float64) # ============================================================================== @@ -267,7 +429,7 @@ julia> mp1 + 5 5.0 ``` """ -mp1 = mpone(Float64) +const global mp1 = mpone(Float64) # ============================================================================== @@ -278,7 +440,7 @@ Max-plus constant "top" equal to +Inf. Equivalent to Scicoslab code: `%top = +Inf` """ -mptop = MP{Float64}(Inf) +const global mptop = MP{Float64}(Inf) # ============================================================================== @@ -295,7 +457,7 @@ not `one()` and as consequence the max-plus identity matrix is not well formed. This const allows to be more algebra compliant by calling `one()` and fixing the fucntion `Matrix{T}(I, m, n)`. """ -const mpI = UniformScaling(one(MP{Float64}).λ) +const global mpI = UniformScaling(one(MP{Float64}).λ) """ mpeye(::Type{T}, n::Int64) @@ -368,7 +530,7 @@ Construct a max-plus one n-by-n matrix. # Examples ```julia-repl -julia> mpones(Float64,2) +julia> mpones(Float64, 2) 2-element Array{MP{Float64},1}: 0.0 0.0 @@ -385,7 +547,7 @@ Construct a max-plus one m-by-n matrix. # Examples ```julia-repl -julia> mpones(Float64,2) +julia> mpones(Float64, 2,2) 2×2 Array{MP{Float64},2}: 0.0 0.0 0.0 0.0 @@ -396,21 +558,58 @@ mpones(::Type{T}, m::Int64, n::Int64) where T = ones(MP{T}, m, n) # ============================================================================== """ - mpsparse(A::Array{T}) + mpsparse(A::Array{T}; preserve=false) Transform a dense matrix to a max-plus sparse matrix. + +# Examples +```julia-repl +julia> mpsparse([4 0; 7 -Inf]) +2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries: + [1, 1] = 4.0 + [2, 1] = 7.0 + [1, 2] = 0.0 +``` """ -mpsparse(A::Array{T}) where T = mpsparse(mparray(A)) -function mpsparse(M::Array{MP{T}}) where T - R = Array{Int64, 1}[]; - C = Array{Int64, 1}[]; - V = Array{MP{T}, 1}[]; +mpsparse(A::Array{T}) where T = mpsparse(MP(A)) +# ; preserve=false) where T = mpsparse(MP(A), preserve) - for i = 1:size(M, 1), j=1:size(M, 2) - M[i,j] != zero(MP{T}) && (R = [R; i]; C = [C; j]; V = [V; M[i,j]]) - end +""" + mpsparse(A::Array{MP{T}}; preserve=false) - sparse(R, C, V) +Transform a dense matrix to a max-plus sparse matrix. By default +Remove -Inf values. To keep them set param preserve to true. + +# Examples +```julia-repl +julia> mpsparse([4 0; 7 -Inf]) +2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries: + [1, 1] = 4.0 + [2, 1] = 7.0 + [1, 2] = 0.0 + +mpsparse([4 0; 7 -Inf], preserve=true) +2×2 SparseMatrixCSC{MP{Float64},Int64} with 3 stored entries: + [1, 1] = 4.0 + [2, 1] = 7.0 + [1, 2] = 0.0 + [2, 2] = -Inf +``` +""" +function mpsparse(M::Array{MP{T}}) where T #; preserve=false) where T + #if preserve + # sparse(M) + #else + R = Array{Int64, 1}[]; + C = Array{Int64, 1}[]; + V = Array{MP{T}, 1}[]; + + for j = 1:size(M, 2), i=1:size(M, 1) + M[i,j] != zero(MP{T}) && (R = [R; i]; C = [C; j]; V = [V; M[i,j]]) + end + + sparse(R, C, V) + #end end # ============================================================================== @@ -447,7 +646,11 @@ julia> array(mpzeros(Float64, 2,2)) [2, 2] = -Inf ``` """ -array(A::SparseMatrixCSC{MP{T},Int64}) where T = map(x -> x.λ, A) +array(A::SparseMatrixCSC{MP{T},U}) where T where U = map(x -> x.λ, A) + +# ============================================================================== + +SparseArrays.sparse(S::SparseMatrixCSC{MP{T},U}) where T where U = map(x -> x.λ, S) # ============================================================================== @@ -465,7 +668,7 @@ julia> full(mpzeros(Float64, 2,5)) -Inf -Inf -Inf -Inf -Inf ``` """ -full(A::SparseMatrixCSC{MP{T},Int64}) where T = Matrix(A) +full(S::SparseMatrixCSC{MP{T},U}) where T where U = Matrix(S) # ============================================================================== @@ -483,12 +686,12 @@ julia> dense(mpzeros(Float64, 2,5)) -Inf -Inf -Inf -Inf -Inf ``` """ -dense(A::SparseMatrixCSC{MP{T},Int64}) where T = Matrix(A) +dense(S::SparseMatrixCSC{MP{T},U}) where T where U = Matrix(S) # ============================================================================== """ - plustimes(a::MP{T}) + plustimes(x::MP{T}) Convert max-plus number to the standard type. Note: the function name comes from Scicoslab. @@ -504,26 +707,27 @@ julia> plustimes([MP(1.0) 2.0; 3.0 4.0]) 3.0 4.0 ``` """ -plustimes(a::MP{T}) where T = a.λ -plustimes(a::Array{MP{T}}) where T = array(a) +plustimes(n::MP{T}) where T = n.λ +plustimes(A::Array{MP{T}}) where T = array(A) +plustimes(S::SparseMatrixCSC{MP{T},U}) where T where U = array(S) # ============================================================================== +# Conversion to min-plus algebra. Convert +Inf and -Inf """ - mparray(A::Array) + minplus(A::Array) -Convert an array to a max-plus array. +Conversion to min-plus algebra. Convert +Inf and -Inf # Examples ```julia-repl -julia> mparray([1.0 2; 3 4]) -2×2 Array{MP{Float64},2}: - 1.0 2.0 - 3.0 4.0 +julia> A = MP([0 3 Inf 1; 1 2 2 -Inf; -Inf Inf 1 0]) +minplus(A); ``` """ -mparray(A::Array) = map(MP, A) -# FIXME missing mparray([1, 2, 3], [1, 2, 3], [0, 2, 0]) +minplus(n::MP{T}) where T = map(x -> (x.λ == mp0) ? mptop : ((x.λ == mptop) ? mp0 : x), n) +minplus(A::Array{MP{T}}) where T = map(x -> minplus(x), A) +minplus(S::SparseMatrixCSC{MP{T},U}) where T where U = map(x -> minplus(x), S) # ============================================================================== @@ -535,12 +739,12 @@ Compute the trace of the matrix (summation of diagonal elements). # Examples ```julia-repl julia> mptrace([MP(1) 2; 3 4]) -MP{Float64}(4.0) +4 ``` """ -mptrace(A::Array{MP{T}}) where T = sum(diag(A)) -mptrace(A::Array{T}) where T = sum(diag(mparray(A))) -mptrace(A::SparseMatrixCSC{MP{T},Int64}) where T = sum(diag(A)) -mptrace(A::SparseMatrixCSC{T,Int64}) where T = sum(diag(A)) +mptrace(A::Array{MP{T}}) where T = isempty(A) ? mp0 : sum(diag(A)) +mptrace(A::Array{T}) where T = isempty(A) ? mp0 : sum(MP(diag(A))) +mptrace(S::SparseMatrixCSC{MP{T}}) where T = isempty(S) ? mp0 : sum(diag(S)) +mptrace(S::SparseMatrixCSC{T}) where T = sum(MP(diag(S))) end # MaxPlus module diff --git a/test/runtests.jl b/test/runtests.jl index b572d6d..8db7af1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,8 @@ -using MaxPlus -using Test -using SparseArrays +using + MaxPlus, Test, Base, LinearAlgebra, SparseArrays # ============================================================================== -# Check Typeof and promotions +# Check type of max-plus and type promotions # ============================================================================== a = MP(1.0) @@ -17,7 +16,7 @@ c = MP{Float64}(1) A = [1 2; 3 4] @test typeof(A) == Array{Int64,2} -@test typeof(mparray(A)) == Array{MP{Int64},2} +@test typeof(MP(A)) == Array{MP{Int64},2} B = [MP(1) MP(2); MP(3) MP(4)] @test typeof(B) == Array{MP{Int64},2} @@ -38,13 +37,23 @@ F = [2 MP(1.0); 3 4] @test typeof(F) == Array{MP{Float64},2} @test typeof(array(F)) == Array{Float64,2} -@test typeof(array(mparray(A))) == typeof(A) -@test typeof(mparray(array(B))) == typeof(B) +@test typeof(array(MP(A))) == typeof(A) +@test typeof(MP(array(B))) == typeof(B) + +# ============================================================================== +# Max-plus constructor +# ============================================================================== + +c = MP(1.0) +@test typeof(MP(c)) == MP{Float64} +d = MP(c) +@test typeof(d) == MP{Float64} # ============================================================================== # Max-plus comparaisons # ============================================================================== +# Max-plus scalars b = MP(3.0) @test (b == b) == true @test (b != b) == (b ≠ b) == false @@ -53,17 +62,25 @@ b = MP(3.0) @test (b > b) == false @test (b < b) == false -# TODO: matrix of bool -# B = [MP(1) MP(2); MP(3) MP(4)] +# Max-plus vs. real comparaison +@test (b == 3.0) == true +@test (b == 4.0) == false +@test (b != 4.0) == true +@test (b < 4.0) == true +@test (b <= 4.0) == true +@test (b > 4.0) == false +@test (b >= 4.0) == false +@test (4.0 > b) == true +@test (4.0 >= b) == true + +# Matrix comparaison +B = [MP(1) MP(2); MP(3) MP(4)] @test (B == B) == true @test (B != B) == (B ≠ B) == false -# @test (B >= B) == (B ≥ B) == true -# @test (B <= B) == (B ≤ B) == true -# @test (B > B) == false -# @test (B < B) == false - -@test array(mparray(A)) == A -@test mparray(array(B)) == B +@test (B .>= B) == (B .≥ B) == [true true; true true] +@test (B .<= B) == (B .≤ B) == [true true; true true] +@test (B .> B) == [false false; false false] +@test (B .< B) == [false false; false false] # ============================================================================== # Max-Plus zero @@ -73,6 +90,7 @@ b = MP(3.0) @test mp0 == zero(MP{Float64}) @test zero(MP{Int64}) == -9223372036854775808 @test zero(MP{Bool}) == false +@test zero(MP(42.0)) == mp0 # ============================================================================== # Max-Plus one @@ -82,6 +100,7 @@ b = MP(3.0) @test mp1 == one(MP{Float64}) @test one(MP{Int64}) == 0 @test one(MP{Bool}) == false +@test one(MP(42.0)) == mp1 # ============================================================================== # Max-plus operations and neutral elements and commutativity @@ -99,8 +118,6 @@ b = MP(3.0) @test a + b == MP(max(5.0, 3.0)) == b + a == MP(max(3.0, 5.0)) == MP(5.0) @test a * b == MP(5.0 + 3.0) == b * a == MP(3.0 + 5.0) == MP(8.0) -@test MP(2)^3 == MP(2 * 3) == MP(6) - # ============================================================================== # Max-plus distributivity of operations # ============================================================================== @@ -116,10 +133,24 @@ c = MP(1.0) @test (a + b) * c == MP(max(5.0, 3.0) + 1.0) == MP(6.0) @test (a * c) + (b * c) == MP(max(5.0 + 1.0, 3.0 + 1.0)) == MP(6.0) +# ============================================================================== +# Max-plus element by element operations +# ============================================================================== + +A=MP([1.0 2.0; 3.0 4.0]) +@test 2.0 .+ A == A .+ 2.0 == [MP(2.0) + MP(1.0) MP(2.0) + MP(2.0); MP(2.0) + MP(3.0) MP(2.0) + MP(4.0)] == MP([2.0 2.0; 3.0 4.0]) +@test 2.0 .* A == A .* 2.0 == [MP(2.0) * MP(1.0) MP(2.0) * MP(2.0); MP(2.0) * MP(3.0) MP(2.0) * MP(4.0)] == MP([3.0 4.0; 5.0 6.0]) + # ============================================================================== # Max-plus other operations # ============================================================================== +@test MP(2)^4 == MP(2 * 4) == MP(8) +@test MP(0)^0 == MP(0 * 0) == MP(0) +@test MP(2)^0 == MP(2 * 0) == MP(0) +@test MP(2)^-3 == MP(2)^(-3) == MP(2 * -3) == MP(-6) +@test MP(2)^0.5 == MP(2 * 0.5) == MP(1.0) + @test b / b == MP(3.0 - 3.0) == MP(0.0) @test b - b == MP(3.0 - 3.0) == MP(0.0) @@ -142,12 +173,10 @@ b = MP(3.0) # Max-plus dense matrix creation # ============================================================================== -b = MP(3.0); C = [b 4.0; 5.0 6.0] +b = MP(3.0); C = [b 4.0; 5.0 6.0]; D = MP([3.0 4.0; 5.0 6.0]) +@test C == D @test typeof(C) == Array{MP{Float64},2} -D = mparray([1.0 2.0; 3.0 4.0]) -@test typeof(D) == Array{MP{Float64},2} - # Max-Plus matrix of max-plus ones O = mpones(Float64, 2) @test typeof(O) == Array{MP{Float64},1} @@ -158,49 +187,93 @@ O = mpones(Float64, 2,5) @test O == [mp1 mp1 mp1 mp1 mp1; mp1 mp1 mp1 mp1 mp1] # Identity matrix +Id = mpeye(Float64, 2) +@test Id == mpeye(Float64, 2) +@test typeof(Id) == Array{MP{Float64},2} +@test Id == [mp1 mp0; mp0 mp1] + Id = mpeye(Float64, 2,5) @test typeof(Id) == Array{MP{Float64},2} @test Id == [mp1 mp0 mp0 mp0 mp0; mp0 mp1 mp0 mp0 mp0] -# TODO -# E = mparray([1, 2, 3], [1, 2, 3], [0, 2, 0]) +# ============================================================================== +# Max-plus / non max-plus conversion +# ============================================================================== + +n = 2.0 +@test plustimes(MP(n)) == n +@test plustimes(mpzeros(Float64,2,2)) == ones(Float64, 2,2) .* mp0 +@test array(MP(A)) == plustimes(MP(A)) == A +@test MP(array(B)) == MP(plustimes(B)) == B + +# ============================================================================== +# Convert values usable for min-plus +# ============================================================================== +@test minplus(MP(1.0)) == MP(1.0) +@test minplus(mp0) == mptop +@test minplus(mptop) == mp0 +@test minplus(MP([Inf 0.0; 7 -Inf])) == MP([-Inf 0.0; 7 Inf]) +#@test minplus(mpsparse([Inf 0.0; 7 -Inf], preserve=true)) == mpsparse([-Inf 0.0; 7 Inf], preserve=true) +A = MP([0 3 Inf 1; 1 2 2 -Inf; -Inf Inf 1 0]) +@test minplus(A) == MP([0 3 -Inf 1; 1 2 2 Inf; Inf -Inf 1 0]) +@test minplus(minplus(A)) == A # ============================================================================== # Max-plus sparse matrix creation and conversion # ============================================================================== +# Construct a sparse max-plus vector +V = MP(sparse([1.0, 0.0, 1.0])) +@test typeof(V) == SparseVector{MP{Float64},Int64} +@test V.nzval == MP([1.0; 1.0]) + +# Construct a sparse max-plus vector of zeros +V = mpzeros(Float64, 2) +@test typeof(V) == SparseVector{MP{Float64},Int64} +@test V.nzval == MP([]) + +# Construct a sparse max-plus matrix +A = MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0])) +@test typeof(A) == SparseMatrixCSC{MP{Float64},Int64} +@test A.nzval == [MP(-Inf), MP(2.0), MP(0.0)] + +# Construct a sparse max-plus matrix +A = MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]), preserve=false) +@test typeof(A) == SparseMatrixCSC{MP{Float64},Int64} +@test A.nzval == [MP(2.0)] + # Basic dense non max-plus matrix to max-plus sparse array A = mpsparse([1.0 2.0; 3.0 4.0]) @test typeof(A) == SparseMatrixCSC{MP{Float64},Int64} @test (A.m == 2) && (A.n == 2) -@test nonzeros(A) == mparray([1.0; 3.0; 2.0; 4.0]) +@test nonzeros(A) == MP([1.0; 3.0; 2.0; 4.0]) -A = mpsparse(mparray([1.0 2.0; 3.0 4.0])) +A = mpsparse(MP([1.0 2.0; 3.0 4.0])) @test typeof(A) == SparseMatrixCSC{MP{Float64},Int64} @test (A.m == 2) && (A.n == 2) -@test nonzeros(A) == mparray([1.0; 3.0; 2.0; 4.0]) +@test nonzeros(A) == MP([1.0; 3.0; 2.0; 4.0]) # 0.0 are not eliminated in the max-plus sparse array B = mpsparse([1.0 2.0; 0.0 4.0]) @test typeof(B) == SparseMatrixCSC{MP{Float64},Int64} @test (B.m == 2) && (B.n == 2) -@test nonzeros(B) == mparray([1.0; 0.0; 2.0; 4.0]) +@test nonzeros(B) == MP([1.0; 0.0; 2.0; 4.0]) -B = mpsparse(mparray([1.0 2.0; 0.0 4.0])) +B = mpsparse(MP([1.0 2.0; 0.0 4.0])) @test typeof(B) == SparseMatrixCSC{MP{Float64},Int64} @test (B.m == 2) && (B.n == 2) -@test nonzeros(B) == mparray([1.0; 0.0; 2.0; 4.0]) +@test nonzeros(B) == MP([1.0; 0.0; 2.0; 4.0]) # -Inf are eliminated in the max-plus sparse array C = mpsparse([1.0 2.0; -Inf 4.0]) @test typeof(C) == SparseMatrixCSC{MP{Float64},Int64} @test (C.m == 2) && (C.n == 2) -@test nonzeros(C) == mparray([1.0; 2.0; 4.0]) +@test nonzeros(C) == MP([1.0; 2.0; 4.0]) -C = mpsparse(mparray([1.0 2.0; -Inf 4.0])) +C = mpsparse(MP([1.0 2.0; -Inf 4.0])) @test typeof(C) == SparseMatrixCSC{MP{Float64},Int64} @test (C.m == 2) && (C.n == 2) -@test nonzeros(C) == mparray([1.0; 2.0; 4.0]) +@test nonzeros(C) == MP([1.0; 2.0; 4.0]) # Max-Plus matrix of zeros is already a max-plus sparse array D = mpzeros(Float64, 3,4) @@ -212,7 +285,7 @@ D = mpzeros(Float64, 3,4) E = mpsparse(mpones(Float64, 3,4)) @test typeof(E) == SparseMatrixCSC{MP{Float64},Int64} @test (E.m == 3) && (E.n == 4) -@test nonzeros(E) == mparray(zeros(Float64, 12)) +@test nonzeros(E) == MP(zeros(Float64, 12)) # Convert non max-plus array to max-plus sparse array F = mpsparse(ones(Float64, 3,4)) @@ -221,10 +294,10 @@ F = mpsparse(ones(Float64, 3,4)) @test nonzeros(F) == ones(12) .+ MP(1.0) # -G = mpsparse(mparray([1.0 2.0 3.0; 1.0 2.0 3.0; 0.0 2.0 0.0])) +G = mpsparse(MP([1.0 2.0 3.0; 1.0 2.0 3.0; 0.0 2.0 0.0])) @test typeof(G) == SparseMatrixCSC{MP{Float64},Int64} @test (G.m == 3) && (G.n == 3) -@test nonzeros(G) == mparray([1.0; 1.0; 0.0; 2.0; 2.0; 2.0; 3.0; 3.0; 0.0]) +@test nonzeros(G) == MP([1.0; 1.0; 0.0; 2.0; 2.0; 2.0; 3.0; 3.0; 0.0]) # max-plus sparse array to max-plus dense array Z = dense(mpzeros(Float64, 2,2)) @@ -245,21 +318,63 @@ Z = array(mpzeros(Float64, 2,2)) # Max-plus operations with matrices # ============================================================================== -A = mparray([4 3; 7 -Inf]) -@test A * A == A^2 == mparray([10.0 7.0; 11.0 10.0]) +A = MP([4 3; 7 -Inf]) +@test A * A == A^2 == MP([10.0 7.0; 11.0 10.0]) C = [MP(3.0) 4.0; 5.0 6.0] -D = mparray([1.0 2; 3 4]) -@test C + D == D + C == mparray([3.0 4.0; 5.0 6.0]) == C +D = MP([1.0 2; 3 4]) +@test C + D == D + C == MP([3.0 4.0; 5.0 6.0]) == C # Trace -@test mptrace(mpeye(Float64, 2,2)) == 0.0 +# TODO shall we forbid trace of non squared matrices ? +@test mptrace([]) == mp0 +@test mptrace(mpeye(Float64, 2,2)) == tr(mpeye(Float64, 2,2)) == 0.0 @test mptrace(mpeye(Float64, 2,5)) == 0.0 -@test mptrace(mpzeros(Float64, 2,2)) == mptrace(full(mpzeros(Float64, 2,2))) == mp0 +@test mptrace(mpzeros(Float64, 2,2)) == mptrace(full(mpzeros(Float64, 2,2))) == tr(mpzeros(Float64, 2,2)) == mp0 @test mptrace(mpzeros(Float64, 2,5)) == mptrace(full(mpzeros(Float64, 2,5))) == mp0 -@test mptrace([1.0 2.0; 3.0 4.0]) == MP(1.0) + MP(4.0) == MP(4.0) +@test mptrace([1.0 2.0; 3.0 4.0]) == tr(MP([1.0 2.0; 3.0 4.0])) == MP(1.0) + MP(4.0) == MP(4.0) +@test mptrace(sparse([1.0 2.0; 3.0 4.0])) == mptrace(mpsparse([1.0 2.0; 3.0 4.0])) == MP(4.0) -# FIXME @test C / C == mparray([0.0 -2.0; 2.0 0.0]) +# ============================================================================== +# +# ============================================================================== + +# FIXME @test C / C == MP([0.0 -2.0; 2.0 0.0]) + +# ============================================================================== +# Max-plus display +# ============================================================================== +using Suppressor + +mp_set_show_mode(0) +result = @capture_out show(stdout, mp0) +@test result == "-Inf" +result = @capture_out show(stdout, mp1) +@test result == "0.0" +result = @capture_out show(stdout, MP(4)) +@test result == "4" + +mp_set_show_mode(1) +result = @capture_out show(stdout, mp0) +@test result == "ε" +result = @capture_out show(stdout, mp1) +@test result == "0.0" +result = @capture_out show(stdout, MP(4)) +@test result == "4" + +mp_set_show_mode(2) +result = @capture_out show(stdout, mp0) +@test result == "ε" +result = @capture_out show(stdout, mp1) +@test result == "e" +result = @capture_out show(stdout, MP(4)) +@test result == "4" + +A = MP([4 0.0; 7 -Inf]) +result = @capture_out show(stdout, A) +@test result == "MP{Float64}[4.0 e; 7.0 ε]" +result = @capture_out LaTeX(stdout, A) +@test result == "\\left[\n\\begin{array}{*{20}c}\n4 & 7 \\\\\ne & \\varepsilon \\\\\n\\end{array}\n\\right]\n" # ============================================================================== #