Don't error if the loglikelihood isn't finite (#127)

* Don't error if the loglikelihood isn't finite

* Fix

* Fix

Author: gdalle

Project.toml

@@ -1,7 +1,7 @@
 name = "HiddenMarkovModels"
 uuid = "84ca31d5-effc-45e0-bfda-5a68cd981f47"
 authors = ["Guillaume Dalle"]
-version = "0.6.1"
+version = "0.6.2"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"

src/inference/forward.jl

@@ -68,11 +68,12 @@ function _forward_digest_observation!(
     current_obs_likelihoods::AbstractVector{<:Real},
     hmm::AbstractHMM,
     obs,
-    control,
+    control;
+    error_if_not_finite::Bool,
 )
     a, b = current_state_marginals, current_obs_likelihoods
 
-    obs_logdensities!(b, hmm, obs, control)
+    obs_logdensities!(b, hmm, obs, control; error_if_not_finite)
     logm = maximum(b)
     b .= exp.(b .- logm)
 
@@ -90,7 +91,8 @@ function _forward!(
     obs_seq::AbstractVector,
     control_seq::AbstractVector,
     seq_ends::AbstractVectorOrNTuple{Int},
-    k::Integer,
+    k::Integer;
+    error_if_not_finite::Bool,
 )
     (; α, B, c, logL) = storage
     t1, t2 = seq_limits(seq_ends, k)
@@ -104,12 +106,14 @@ function _forward!(
             αₜ₋₁ = view(α, :, t - 1)
             predict_next_state!(αₜ, hmm, αₜ₋₁, control_seq[t - 1])
         end
-        cₜ, logLₜ = _forward_digest_observation!(αₜ, Bₜ, hmm, obs_seq[t], control_seq[t])
+        cₜ, logLₜ = _forward_digest_observation!(
+            αₜ, Bₜ, hmm, obs_seq[t], control_seq[t]; error_if_not_finite
+        )
         c[t] = cₜ
         logL[k] += logLₜ
     end
 
-    @argcheck isfinite(logL[k])
+    error_if_not_finite && @argcheck isfinite(logL[k])
     return nothing
 end
 
@@ -122,14 +126,15 @@ function forward!(
     obs_seq::AbstractVector,
     control_seq::AbstractVector;
     seq_ends::AbstractVectorOrNTuple{Int},
+    error_if_not_finite::Bool=true,
 )
     if seq_ends isa NTuple{1}
         for k in eachindex(seq_ends)
-            _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k)
+            _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
         end
     else
         @threads for k in eachindex(seq_ends)
-            _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k)
+            _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
         end
     end
     return nothing
@@ -147,8 +152,9 @@ function forward(
     obs_seq::AbstractVector,
     control_seq::AbstractVector=Fill(nothing, length(obs_seq));
     seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),),
+    error_if_not_finite::Bool=true,
 )
     storage = initialize_forward(hmm, obs_seq, control_seq; seq_ends)
-    forward!(storage, hmm, obs_seq, control_seq; seq_ends)
+    forward!(storage, hmm, obs_seq, control_seq; seq_ends, error_if_not_finite)
     return storage.α, storage.logL
 end

src/inference/forward_backward.jl

@@ -42,7 +42,7 @@ function _forward_backward!(
     t1, t2 = seq_limits(seq_ends, k)
 
     # Forward (fill B, α, c and logL)
-    _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k)
+    _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite=true)
 
     # Backward
     β[:, t2] .= c[t2]

src/inference/logdensity.jl

@@ -9,7 +9,7 @@ function DensityInterface.logdensityof(
     control_seq::AbstractVector=Fill(nothing, length(obs_seq));
     seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),),
 )
-    _, logL = forward(hmm, obs_seq, control_seq; seq_ends)
+    _, logL = forward(hmm, obs_seq, control_seq; seq_ends, error_if_not_finite=false)
     return sum(logL)
 end
 

src/types/abstract_hmm.jl

@@ -128,13 +128,13 @@ StatsAPI.fit!
 ## Fill logdensities
 
 function obs_logdensities!(
-    logb::AbstractVector{T}, hmm::AbstractHMM, obs, control
+    logb::AbstractVector{T}, hmm::AbstractHMM, obs, control; error_if_not_finite::Bool=true
 ) where {T}
     dists = obs_distributions(hmm, control)
     @simd for i in eachindex(logb, dists)
         logb[i] = logdensityof(dists[i], obs)
     end
-    @argcheck maximum(logb) < typemax(T)
+    error_if_not_finite && @argcheck maximum(logb) < typemax(T)
     return nothing
 end
 

test/misc.jl

@@ -0,0 +1,13 @@
+using HiddenMarkovModels
+using HiddenMarkovModels: rand_prob_vec, rand_trans_mat
+using Distributions
+using Test
+
+@testset "Allow NaN density" begin
+    init = rand_prob_vec(2)
+    trans = rand_trans_mat(2)
+    dists = [Normal(i, Inf) for i in 1:2]
+    hmm = HMM(init, trans, dists)
+    obs_seq = rand(5)
+    @test isnan(logdensityof(hmm, obs_seq))
+end

test/runtests.jl

@@ -53,4 +53,8 @@ Pkg.develop(; path=joinpath(dirname(@__DIR__), "libs", "HMMTest"))
     @testset verbose = true "Correctness - $TEST_SUITE" begin
         include("correctness.jl")
     end
+
+    @testset verbose = true "Miscellaneous" begin
+        include("misc.jl")
+    end
 end