JuliaArrays / StaticArrays.jl

Statically sized arrays for Julia
Other
764 stars 148 forks source link

NaN and Inf issues in exp #875

Open andreasnoack opened 3 years ago

andreasnoack commented 3 years ago

The test in https://github.com/JuliaArrays/StaticArrays.jl/blob/59f92e0ca7ac391a850a6e7a2ce1eb53aa237fc4/src/expm.jl#L21-L32 don't cover the NaN case so

julia> A = @SMatrix([1.0 0.1; NaN 1.0])
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.1
 NaN    1.0

julia> exp(A)
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
   2.71828  0.271828
 NaN        2.71828

which should either throw or give a NaN. Also, https://github.com/JuliaArrays/StaticArrays.jl/blob/59f92e0ca7ac391a850a6e7a2ce1eb53aa237fc4/src/expm.jl#L96-L99 can be reached for matrices with Infs which will error in an uninformative way

julia> A = @SMatrix([1.0 0.1 0.0; Inf 1.0 0.1; 0.1 0.0 1.0])
3×3 SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
  1.0  0.1  0.0
 Inf   1.0  0.1
  0.1  0.0  1.0

julia> exp(A)
ERROR: InexactError: trunc(Int64, Inf)
Stacktrace:
 [1] trunc at ./float.jl:703 [inlined]
 [2] ceil at ./float.jl:365 [inlined]
 [3] _exp(::Size{(3, 3)}, ::SArray{Tuple{3,3},Float64,2,9}) at /Users/vagrant/.julia/packages/StaticArrays/LJQEe/src/expm.jl:97
 [4] exp(::SArray{Tuple{3,3},Float64,2,9}) at /Users/vagrant/.julia/packages/StaticArrays/LJQEe/src/expm.jl:1
 [5] top-level scope at REPL[70]:1
 [6] eval(::Module, ::Any) at /Applications/Pumas-1.1.app/Contents/Resources/julia/Contents/Resources/julia/lib/julia/sys.dylib:?
 [7] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/vagrant/worker/juliapro-staging-osx1011-0_6/build/tmp_julia/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/REPL/src/REPL.jl:134
 [8] repl_backend_loop(::REPL.REPLBackend) at /Users/vagrant/worker/juliapro-staging-osx1011-0_6/build/tmp_julia/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/REPL/src/REPL.jl:195
 [9] start_repl_backend(::REPL.REPLBackend, ::Any) at /Users/vagrant/worker/juliapro-staging-osx1011-0_6/build/tmp_julia/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/REPL/src/REPL.jl:180
 [10] run_repl(::REPL.AbstractREPL, ::Any; backend_on_current_task::Bool) at /Users/vagrant/worker/juliapro-staging-osx1011-0_6/build/tmp_julia/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/REPL/src/REPL.jl:292
 [11] run_repl(::REPL.AbstractREPL, ::Any) at /Users/vagrant/worker/juliapro-staging-osx1011-0_6/build/tmp_julia/Julia-1.5.app/Contents/Resources/julia/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
 [12] (::Base.var"#807#809"{Bool,Bool,Bool,Bool})(::Module) at ./client.jl:399
 [13] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at /Applications/Pumas-1.1.app/Contents/Resources/julia/Contents/Resources/julia/lib/julia/sys.dylib:?
 [14] exec_options(::Base.JLOptions) at /Applications/Pumas-1.1.app/Contents/Resources/julia/Contents/Resources/julia/lib/julia/sys.dylib:?
 [15] _start() at /Applications/Pumas-1.1.app/Contents/Resources/julia/Contents/Resources/julia/lib/julia/sys.dylib:?
andreasnoack commented 2 years ago

Adding one more example here

julia> A = @SMatrix([-1.0 0.0; 1.0 -0.1875])
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 -1.0   0.0
  1.0  -0.1875

julia> exp(2016*A)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 NaN  NaN
 NaN  NaN

julia> exp(2016*Matrix(A))
2×2 Matrix{Float64}:
 0.0           0.0
 8.45011e-165  6.86572e-165