diff --git a/test/arhmm_testing_plutonotebook.jl b/test/arhmm_testing_plutonotebook.jl deleted file mode 100644 index 4e4ab82..0000000 --- a/test/arhmm_testing_plutonotebook.jl +++ /dev/null @@ -1,1638 +0,0 @@ -### A Pluto.jl notebook ### -# v0.20.3 - -using Markdown -using InteractiveUtils - -# ╔═╡ 455866a4-a5ef-11ef-04ff-f347c85da0d8 -begin - using ArgCheck: @argcheck - using Base: RefValue - using Base.Threads: @threads - using ChainRulesCore: ChainRulesCore, NoTangent, RuleConfig, rrule_via_ad - using DensityInterface: DensityInterface, DensityKind, HasDensity, NoDensity, logdensityof - using DocStringExtensions - using FillArrays: Fill - using LinearAlgebra: Transpose, axpy!, dot, ldiv!, lmul!, mul!, parent - using Random: Random, AbstractRNG, default_rng - using SparseArrays: AbstractSparseArray, SparseMatrixCSC, nonzeros, nnz, nzrange, rowvals - using StatsAPI: StatsAPI, fit, fit! - using StatsFuns: log2π -end - -# ╔═╡ b8b4c6ec-a8f3-4609-a0c0-15652c0d15be -begin - using Distributions - using PlutoUI - TableOfContents() -end - -# ╔═╡ a5c89117-9642-4c28-bf5e-7c93184934b4 -md""" -## 1. abstract_hmm.jl -""" - -# ╔═╡ 5d513c86-0595-4fee-8312-9927ea731fde -md""" -## 2. linalg.jl -""" - -# ╔═╡ b6abd844-345b-4c36-b6c0-699fa5c45fd9 -md""" -## 3. limits.jl -""" - -# ╔═╡ ebfa27fc-bdb0-4aaa-8930-85025df33b7c -md""" -## 4. forward.jl -""" - -# ╔═╡ f2599655-855a-4c41-a6ed-c22c16d792a1 -md""" -## 5. logdensity.jl -""" - -# ╔═╡ e35f3103-c5ec-4a6f-bfc0-76326f8eed54 -md""" -## 6. forward_backward.jl.jl -""" - -# ╔═╡ 71193717-7f4e-4db6-9669-149e7d6405af -md""" -## 7. viterbi.jl -""" - -# ╔═╡ ffb943ac-b403-4f71-a19b-0875748674ed -md""" -## 8. Test -""" - -# ╔═╡ 3a00c710-e6e6-4b5e-875e-7bbcb173d8a3 -# Supposing control Uₜ, |Uₜ| = 2, state Xₜ, |Xₜ| = 2, and observation Yₜ, |Yₜ| = 3 - -# ╔═╡ 8d99e582-fed1-42bd-8196-5f73cffe5a05 -# Discrete ARHMM without control - -# ╔═╡ c3ed40cc-0921-411e-aa8b-547d35fa1842 -uᵢ = [0.5, 0.5] - -# ╔═╡ c24abed8-65bc-4bb8-9444-64377c90be21 -tᵢⱼ = [0.3 0.7; 0.4 0.6] - -# ╔═╡ 4b6ecf0d-6330-4fc5-b31a-f49bbc03d397 -eʳᵢⱼ = [[0.25 0.5 0.25; 0.4 0.4 0.2], [0.7 0.2 0.1; 0.3 0.3 0.4], [0.4 0.2 0.4; 0.1 0.1 0.8]] - -# ╔═╡ 1c825864-2332-4e99-8d3a-4583f4fb56e2 -pᵢⱼ = [0.4 0.4 0.2; 0.3 0.3 0.4] - -# ╔═╡ 3b669bd8-21e1-42d5-b2d6-45340c57e1e2 -# Discrete ARHMM with control - -# ╔═╡ dfa9f848-b928-4299-8347-2ac492958df7 -uᵏᵢ = [[0.1, 0.9], [0.5, 0.5]] - -# ╔═╡ b454b317-29a9-41ec-82a7-9c48a88f2576 -tᵏᵢⱼ = [[0.3 0.7; 0.4 0.6], [0.75 0.25; 0.5 0.5]] - -# ╔═╡ b1377b4b-ce15-4eb8-83cb-740a58673ac6 -eᵏʳᵢⱼ = [[[0.25 0.5 0.25; 0.4 0.4 0.2], [0.7 0.2 0.1; 0.3 0.3 0.4], [0.4 0.2 0.4; 0.1 0.1 0.8]], [[0.3 0.2 0.5; 0.2 0.4 0.4], [0.5 0.25 0.25; 0.5 0.4 0.1], [0.25 0.25 0.5; 0.2 0.6 0.2]]] - -# ╔═╡ e7acfa37-371b-437c-aa78-47e37500d414 -begin - nstates = 2 - ncontrol = 2 - nobservations = 3 - N = 20 - true_control_seq = rand([1, 2], N) - true_state_seq = Int64[] - push!(true_state_seq, rand(Categorical(uᵏᵢ[true_control_seq[1]]))) - for i ∈ 2:N - push!(true_state_seq, rand(Categorical(tᵏᵢⱼ[true_control_seq[i]][true_state_seq[i-1],:]))) - end - true_obs_seq = Int64[] - push!(true_obs_seq, rand(Categorical(eᵏʳᵢⱼ[true_control_seq[1]][rand([1,2,3])][true_state_seq[1],:]))) - for i ∈ 2:N - push!(true_obs_seq, rand(Categorical(eᵏʳᵢⱼ[true_control_seq[i]][true_obs_seq[i-1]][true_state_seq[i],:]))) - end -end - -# ╔═╡ 4d76c441-429f-4500-9c26-6247abad17ed -begin - - const AbstractVectorOrNTuple{T} = Union{AbstractVector{T},NTuple{N,T}} where {N} - - sum_to_one!(x) = ldiv!(sum(x), x) - - mynonzeros(x::AbstractArray) = x - mynonzeros(x::AbstractSparseArray) = nonzeros(x) - - mynnz(x::AbstractArray) = length(mynonzeros(x)) - - elementwise_log(x::AbstractArray) = log.(x) - - function elementwise_log(A::SparseMatrixCSC) - return SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, log.(A.nzval)) - end - - """ - mul_rows_cols!(B, l, A, r) - - Perform the in-place operation `B .= l .* A .* transpose(r)`. - """ - function mul_rows_cols!( - B::AbstractMatrix, l::AbstractVector, A::AbstractMatrix, r::AbstractVector - ) - B .= l .* A .* transpose(r) - return B - end - - function mul_rows_cols!( - B::SparseMatrixCSC, l::AbstractVector, A::SparseMatrixCSC, r::AbstractVector - ) - @argcheck axes(A, 1) == eachindex(r) - @argcheck axes(A, 2) == eachindex(l) - @argcheck size(A) == size(B) - @argcheck nnz(B) == nnz(A) - Brv = rowvals(B) - Bnz = nonzeros(B) - Anz = nonzeros(A) - @simd for j in axes(B, 2) - @argcheck nzrange(B, j) == nzrange(A, j) - @simd for k in nzrange(B, j) - i = Brv[k] - Bnz[k] = l[i] * Anz[k] * r[j] - end - end - return B - end - - """ - argmaxplus_transmul!(y, ind, A, x) - - Perform the in-place multiplication `transpose(A) * x` _in the sense of max-plus algebra_, store the result in `y`, and store the index of the maximum for each component of `y` in `ind`. - """ - function argmaxplus_transmul!( - y::AbstractVector{R}, - ind::AbstractVector{<:Integer}, - A::AbstractMatrix, - x::AbstractVector, - ) where {R} - @argcheck axes(A, 1) == eachindex(x) - @argcheck axes(A, 2) == eachindex(y) - fill!(y, typemin(R)) - fill!(ind, 0) - @simd for j in axes(A, 2) - @simd for i in axes(A, 1) - z = A[i, j] + x[i] - if z > y[j] - y[j] = z - ind[j] = i - end - end - end - return y - end - - function argmaxplus_transmul!( - y::AbstractVector{R}, - ind::AbstractVector{<:Integer}, - A::SparseMatrixCSC, - x::AbstractVector, - ) where {R} - @argcheck axes(A, 1) == eachindex(x) - @argcheck axes(A, 2) == eachindex(y) - Anz = nonzeros(A) - Arv = rowvals(A) - fill!(y, typemin(R)) - fill!(ind, 0) - @simd for j in axes(A, 2) - @simd for k in nzrange(A, j) - i = Arv[k] - z = Anz[k] + x[i] - if z > y[j] - y[j] = z - ind[j] = i - end - end - end - return y - end - -end - -# ╔═╡ 6444a4ab-5992-4a39-a2e9-8301f9a6ba91 -begin - ############################################################################ - # 1. TYPE # - ############################################################################ - """ - AbstractHMM - - Abstract supertype for an HMM amenable to simulation, inference and learning. - - # Interface - - To create your own subtype of `AbstractHMM`, you need to implement the following methods: - - - [`initialization`](@ref) - - [`transition_matrix`](@ref) - - [`obs_distributions`](@ref) - - [`fit!`](@ref) (for learning) - - # Applicable functions - - Any `AbstractHMM` which satisfies the interface can be given to the following functions: - - - [`rand`](@ref) - - [`logdensityof`](@ref) - - [`forward`](@ref) - - [`viterbi`](@ref) - - [`forward_backward`](@ref) - - [`baum_welch`](@ref) (if `[fit!](@ref)` is implemented) - """ - abstract type AbstractHMM{T} end - - @inline DensityInterface.DensityKind(::AbstractHMM) = HasDensity() - - - ############################################################################ - # 2. INTERFACE # - ############################################################################ - - - #------------------------------ 2.1. length -------------------------------# - """ - length(hmm) - - Return the number of states of `hmm`. - """ - Base.length(hmm::AbstractHMM) = length(initialization(hmm)) - - Base.length(hmm::AbstractHMM, control) = length(initialization(hmm, control)) - - Base.length(hmm::AbstractHMM, ::Nothing) = length(initialization(hmm)) - - - #------------------------------ 2.2. eltype -------------------------------# - """ - eltype(hmm, obs, control) - - Return a type that can accommodate forward-backward computations for `hmm` on observations similar to `obs`. - - It is typically a promotion between the element type of the initialization, the element type of the transition matrix, and the type of an observation logdensity evaluated at `obs`. - """ - function Base.eltype(hmm::AbstractHMM, obs, control) - init_type = eltype(initialization(hmm, control)) - trans_type = eltype(transition_matrix(hmm, control)) - dist = obs_distributions(hmm, control, obs)[1] - logdensity_type = typeof(logdensityof(dist, obs)) - return promote_type(init_type, trans_type, logdensity_type) - end - - - #--------------------------- 2.3. initialization --------------------------# - """ - initialization(hmm) - - Return the vector of initial state probabilities for `hmm`. - """ - function initialization end - - initialization(hmm::AbstractHMM) = hmm.init - - initialization(hmm::AbstractHMM, control) = hmm.init - - initialization(hmm::AbstractHMM, ::Nothing) = hmm.init - - - #------------------------ 2.4. log_initialization -------------------------# - """ - log_initialization(hmm) - - Return the vector of initial state log-probabilities for `hmm`. - - Falls back on `initialization`. - """ - log_initialization(hmm::AbstractHMM) = elementwise_log(initialization(hmm)) - - log_initialization(hmm::AbstractHMM, control) = elementwise_log(initialization(hmm, control)) - - log_initialization(hmm::AbstractHMM, ::Nothing) = elementwise_log(initialization(hmm)) - - - #------------------------- 2.5. transition_matrix -------------------------# - """ - transition_matrix(hmm) - transition_matrix(hmm, control) - - Return the matrix of state transition probabilities for `hmm` (possibly when `control` is applied). - - !!! note - When processing sequences, the control at time `t` influences the transition from time `t` to `t+1` (and not from time `t-1` to `t`). - """ - function transition_matrix end - - transition_matrix(hmm::AbstractHMM) = hmm.trans - - transition_matrix(hmm::AbstractHMM, control) = transition_matrix(hmm) - - transition_matrix(hmm::AbstractHMM, ::Nothing) = transition_matrix(hmm) - - - #----------------------- 2.6. log_transition_matrix -----------------------# - """ - log_transition_matrix(hmm) - log_transition_matrix(hmm, control) - - Return the matrix of state transition log-probabilities for `hmm` (possibly when `control` is applied). - - Falls back on `transition_matrix`. - - !!! note - When processing sequences, the control at time `t` influences the transition from time `t` to `t+1` (and not from time `t-1` to `t`). - """ - log_transition_matrix(hmm::AbstractHMM) = elementwise_log(transition_matrix(hmm)) - - log_transition_matrix(hmm::AbstractHMM, control) = elementwise_log(transition_matrix(hmm, control)) - - log_transition_matrix(hmm::AbstractHMM, ::Nothing) = elementwise_log(transition_matrix(hmm)) - - - - #------------------------- 2.7. obs_distributions -------------------------# - """ - obs_distributions(hmm) - obs_distributions(hmm, control) - - Return a vector of observation distributions, one for each state of `hmm` (possibly when `control` is applied). - - These distribution objects should implement - - - `Random.rand(rng, dist)` for sampling - - `DensityInterface.logdensityof(dist, obs)` for inference - - `StatsAPI.fit!(dist, obs_seq, weight_seq)` for learning - """ - function obs_distributions end - - obs_distributions(hmm::AbstractHMM) = [hmm.dists[i] for i ∈ 1:length(hmm)] - - obs_distributions(hmm::AbstractHMM, control, prev_obs) = obs_distributions(hmm, control, prev_obs) - - ### Fallback when it is not autoregressive and there is no control - obs_distributions(hmm::AbstractHMM, ::Nothing) = obs_distributions(hmm) - - ### Fallback when it is autoregressive, but previous observation is Missing - obs_distributions(hmm::AbstractHMM, control, ::Missing) = obs_distributions(hmm, control) - - ### Fallback when it is autoregressive and there is no control, but observation is missing - obs_distributions(hmm::AbstractHMM, ::Nothing, ::Missing) = obs_distributions(hmm) - - - #---------------------------- 2.8. previous_obs ---------------------------# - previous_obs(::AbstractHMM{false}, obs_seq::AbstractVector, t::Integer) = nothing - - previous_obs(::AbstractHMM{true}, obs_seq::AbstractVector, t::Integer) = obs_seq[t - 1] - - - #--------------------------- 2.9. StatsAPI.fit! ---------------------------# - """ - StatsAPI.fit!( - hmm, fb_storage::ForwardBackwardStorage, - obs_seq, [control_seq]; seq_ends, - ) - - Update `hmm` in-place based on information generated during forward-backward. - - This function is allowed to reuse `fb_storage` as a scratch space, so its contents should not be trusted afterwards. - """ - StatsAPI.fit! - - - #------------------------- 2.10. obs_logdensities! ------------------------# - function obs_logdensities!( - logb::AbstractVector{T}, hmm::AbstractHMM, obs, control, prev_obs - ) where {T} - dists = obs_distributions(hmm, control, prev_obs) - @simd for i in eachindex(logb, dists) - logb[i] = logdensityof(dists[i], obs) - end - @argcheck maximum(logb) < typemax(T) - return nothing - end - - - ############################################################################ - # 3. SAMPLING # - ############################################################################ - - # <------------- Didn't touch it yet! - """ - rand([rng,] hmm, T) - rand([rng,] hmm, control_seq) - - Simulate `hmm` for `T` time steps, or when the sequence `control_seq` is applied. - - Return a named tuple `(; state_seq, obs_seq)`. - """ - function Random.rand(rng::AbstractRNG, hmm::AbstractHMM, control_seq::AbstractVector) - T = length(control_seq) - dummy_log_probas = fill(-Inf, length(hmm)) - - init = initialization(hmm, control) - state_seq = Vector{Int}(undef, T) - state1 = rand(rng, LightCategorical(init, dummy_log_probas)) - state_seq[1] = state1 - - @views for t in 1:(T - 1) - trans = transition_matrix(hmm, control_seq[t]) - state_seq[t + 1] = rand( - rng, LightCategorical(trans[state_seq[t], :], dummy_log_probas) - ) - end - - dists1 = obs_distributions(hmm, control_seq[1], missing) - obs1 = rand(rng, dists1[state1]) - obs_seq = Vector{typeof(obs1)}(undef, T) - obs_seq[1] = obs1 - - for t in 2:T - dists = obs_distributions(hmm, control_seq[t], previous_obs(hmm, obs_seq, t)) - obs_seq[t] = rand(rng, dists[state_seq[t]]) - end - return (; state_seq=state_seq, obs_seq=obs_seq) - end - - function Random.rand(hmm::AbstractHMM, control_seq::AbstractVector) - return rand(default_rng(), hmm, control_seq) - end - - function Random.rand(rng::AbstractRNG, hmm::AbstractHMM, T::Integer) - return rand(rng, hmm, Fill(nothing, T)) - end - - function Random.rand(hmm::AbstractHMM, T::Integer) - return rand(hmm, Fill(nothing, T)) - end - - - ############################################################################ - # 4. PRIOR # - ############################################################################ - """ - logdensityof(hmm) - - Return the prior loglikelihood associated with the parameters of `hmm`. - """ - DensityInterface.logdensityof(hmm::AbstractHMM) = false - - - ############################################################################ - # 5. ARHMM EXAMPLE # - ############################################################################ - ## Test scruct for Discrete ARHMM with control - struct DiscreteCARHMM{T<:Number} <: AbstractHMM{true} - # Initial distribution P(X_{1}|U_{1}), one vector for each control - init::Vector{Vector{T}} - # Transition matrix P(X_{t}|X_{t-1}, U_{t}), one matrix for each control - trans::Vector{Matrix{T}} - # Emission matriz P(Y_{t}|X_{t}, U_{t}), one matriz for each control and each possible observation - dists::Vector{Vector{Matrix{T}}} - # Prior Distribution for P(Y_{1}|X_{1}, U_{1}), one matriz for each control - prior::Vector{Matrix{T}} - end - - initialization(hmm::DiscreteCARHMM, control) = hmm.init[control] - - transition_matrix(hmm::DiscreteCARHMM, control) = hmm.trans[control] - - obs_distributions(hmm::DiscreteCARHMM, control, prev_obs) = [Categorical(hmm.dists[control][prev_obs][i,:]) for i in 1:length(hmm, control)] - - obs_distributions(hmm::DiscreteCARHMM, control, ::Missing) = [Categorical(hmm.prior[control][i,:]) for i in 1:length(hmm, control)] - - - ## Test scruct for Discrete ARHMM without control - struct DiscreteARHMM{T<:Number} <: AbstractHMM{true} - # Initial distribution P(X_{1}) - init::Vector{T} - # Transition matrix P(X_{t}|X_{t-1}) - trans::Matrix{T} - # Emission matriz P(Y_{t}|X_{t}) - dists::Vector{Matrix{T}} - # Prior Distribution for P(Y_{1}|X_{1}) - prior::Matrix{T} - end - - obs_distributions(hmm::DiscreteARHMM, ::Nothing, prev_obs) = [Categorical(hmm.dists[prev_obs][i,:]) for i in 1:length(hmm)] - - obs_distributions(hmm::DiscreteARHMM, ::Nothing, ::Missing) = [Categorical(hmm.prior[i,:]) for i in 1:length(hmm)] -end - -# ╔═╡ c90ce230-7640-4337-b6ae-e9369564e8e7 -begin - function predict_next_state!( - next_state_marginals::AbstractVector{<:Real}, - hmm::AbstractHMM, - current_state_marginals::AbstractVector{<:Real}, - control=nothing, - ) - trans = transition_matrix(hmm, control) - mul!(next_state_marginals, transpose(trans), current_state_marginals) - return next_state_marginals - end - - function predict_previous_state!( - previous_state_marginals::AbstractVector{<:Real}, - hmm::AbstractHMM, - current_state_marginals::AbstractVector{<:Real}, - control=nothing, - ) - trans = transition_matrix(hmm, control) - mul!(previous_state_marginals, trans, current_state_marginals) - return previous_state_marginals - end -end - -# ╔═╡ 5e9a1a49-0356-48ec-ad85-ed5c49653d2d -hmm₀ = DiscreteARHMM(uᵢ, tᵢⱼ, eʳᵢⱼ, pᵢⱼ) - -# ╔═╡ dc54a9a0-4fd6-4f64-9067-775d7674ebc1 -""" -$(SIGNATURES) - -Return a tuple `(t1, t2)` giving the begin and end indices of subsequence `k` within a set of sequences ending at `seq_ends`. -""" -function seq_limits(seq_ends::AbstractVectorOrNTuple{Int}, k::Integer) - if k == 1 - return 1, seq_ends[k] - else - return seq_ends[k - 1] + 1, seq_ends[k] - end -end - -# ╔═╡ f4fe5934-b4a3-412c-89e2-688bdffd7dd2 -begin - struct ForwardStorage{R} - "posterior last state marginals `α[i] = ℙ(X[T]=i | Y[1:T])`" - α::Matrix{R} - "one loglikelihood per observation sequence" - logL::Vector{R} - B::Matrix{R} - c::Vector{R} - end - - struct ForwardBackwardStorage{R,M<:AbstractMatrix{R}} - "posterior state marginals `γ[i,t] = ℙ(X[t]=i | Y[1:T])`" - γ::Matrix{R} - "posterior transition marginals `ξ[t][i,j] = ℙ(X[t]=i, X[t+1]=j | Y[1:T])`" - ξ::Vector{M} - "one loglikelihood per observation sequence" - logL::Vector{R} - B::Matrix{R} - α::Matrix{R} - c::Vector{R} - β::Matrix{R} - Bβ::Matrix{R} - end - - Base.eltype(::ForwardBackwardStorage{R}) where {R} = R - - const ForwardOrForwardBackwardStorage{R} = Union{ - ForwardStorage{R},ForwardBackwardStorage{R} - } - - """ - $(SIGNATURES) - """ - function initialize_forward( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - ) - N, T, K = length(hmm, control_seq[1]), length(obs_seq), length(seq_ends) - R = eltype(hmm, obs_seq[1], control_seq[1]) - α = Matrix{R}(undef, N, T) - logL = Vector{R}(undef, K) - B = Matrix{R}(undef, N, T) - c = Vector{R}(undef, T) - return ForwardStorage(α, logL, B, c) - end - - function _forward_digest_observation!( - current_state_marginals::AbstractVector{<:Real}, - current_obs_likelihoods::AbstractVector{<:Real}, - hmm::AbstractHMM, - obs, - control, - prev_obs, - ) - a, b = current_state_marginals, current_obs_likelihoods - - obs_logdensities!(b, hmm, obs, control, prev_obs) - logm = maximum(b) - b .= exp.(b .- logm) - - a .*= b - c = inv(sum(a)) - lmul!(c, a) - - logL = -log(c) + logm - return c, logL - end - - function _forward!( - storage::ForwardOrForwardBackwardStorage, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector, - seq_ends::AbstractVectorOrNTuple{Int}, - k::Integer, - ) - (; α, B, c, logL) = storage - t1, t2 = seq_limits(seq_ends, k) - logL[k] = zero(eltype(logL)) - for t in t1:t2 - αₜ = view(α, :, t) - Bₜ = view(B, :, t) - if t == t1 - copyto!(αₜ, initialization(hmm, control_seq[t])) - else - αₜ₋₁ = view(α, :, t - 1) - predict_next_state!(αₜ, hmm, αₜ₋₁, control_seq[t - 1]) - end - prev_obs = t == t1 ? missing : previous_obs(hmm, obs_seq, t) - cₜ, logLₜ = _forward_digest_observation!( - αₜ, Bₜ, hmm, obs_seq[t], control_seq[t], prev_obs - ) - c[t] = cₜ - logL[k] += logLₜ - end - - @argcheck isfinite(logL[k]) - return nothing - end - - """ - $(SIGNATURES) - """ - function forward!( - storage::ForwardOrForwardBackwardStorage, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - ) - if seq_ends isa NTuple{1} - for k in eachindex(seq_ends) - _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k) - end - else - @threads for k in eachindex(seq_ends) - _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k) - end - end - return nothing - end - - """ - $(SIGNATURES) - - Apply the forward algorithm to infer the current state after sequence `obs_seq` for `hmm`. - - Return a tuple `(storage.α, storage.logL)` where `storage` is of type [`ForwardStorage`](@ref). - """ - function forward( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector=Fill(nothing, length(obs_seq)); - seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),), - ) - storage = initialize_forward(hmm, obs_seq, control_seq; seq_ends) - forward!(storage, hmm, obs_seq, control_seq; seq_ends) - return storage.α, storage.logL - end - - nothing -end - -# ╔═╡ 4bcf1079-13c6-418e-a37e-ce71e5cc74f1 -begin - """ - $(SIGNATURES) - - Run the forward algorithm to compute the loglikelihood of `obs_seq` for `hmm`, integrating over all possible state sequences. - """ - function DensityInterface.logdensityof( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector=Fill(nothing, length(obs_seq)); - seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),), - ) - _, logL = forward(hmm, obs_seq, control_seq; seq_ends) - return sum(logL) - end - - """ - $(SIGNATURES) - - Run the forward algorithm to compute the the joint loglikelihood of `obs_seq` and `state_seq` for `hmm`. - """ - function joint_logdensityof( - hmm::AbstractHMM, - obs_seq::AbstractVector, - state_seq::AbstractVector, - control_seq::AbstractVector=Fill(nothing, length(obs_seq)); - seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),), - ) - R = eltype(hmm, obs_seq[1], control_seq[1]) - logL = zero(R) - for k in eachindex(seq_ends) - t1, t2 = seq_limits(seq_ends, k) - # Initialization - init = initialization(hmm, control_seq[t1]) - logL += log(init[state_seq[t1]]) - # Transitions - for t in t1:(t2 - 1) - trans = transition_matrix(hmm, control_seq[t]) - logL += log(trans[state_seq[t], state_seq[t + 1]]) - end - # Observations - for t in t1:t2 - dists = obs_distributions(hmm, control_seq[t], previous_obs(hmm, obs_seq, t)) - logL += logdensityof(dists[state_seq[t]], obs_seq[t]) - end - end - return logL - end - - nothing -end - -# ╔═╡ 0a903f68-79dc-4a45-b0fe-e5e17eff3b6f -begin - """ - $(SIGNATURES) - """ - function initialize_forward_backward( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - transition_marginals=true, - ) - N, T, K = length(hmm, control_seq[1]), length(obs_seq), length(seq_ends) - R = eltype(hmm, obs_seq[1], control_seq[1]) - trans = transition_matrix(hmm, control_seq[1]) - M = typeof(similar(trans, R)) - - γ = Matrix{R}(undef, N, T) - ξ = Vector{M}(undef, T) - if transition_marginals - for t in 1:T - ξ[t] = similar(transition_matrix(hmm, control_seq[t]), R) - end - end - logL = Vector{R}(undef, K) - B = Matrix{R}(undef, N, T) - α = Matrix{R}(undef, N, T) - c = Vector{R}(undef, T) - β = Matrix{R}(undef, N, T) - Bβ = Matrix{R}(undef, N, T) - return ForwardBackwardStorage{R,M}(γ, ξ, logL, B, α, c, β, Bβ) - end - - function _forward_backward!( - storage::ForwardBackwardStorage{R}, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector, - seq_ends::AbstractVectorOrNTuple{Int}, - k::Integer; - transition_marginals::Bool=true, - ) where {R} - (; α, β, c, γ, ξ, B, Bβ) = storage - t1, t2 = seq_limits(seq_ends, k) - - # Forward (fill B, α, c and logL) - _forward!(storage, hmm, obs_seq, control_seq, seq_ends, k) - - # Backward - β[:, t2] .= c[t2] - for t in (t2 - 1):-1:t1 - Bβ[:, t + 1] .= view(B, :, t + 1) .* view(β, :, t + 1) - βₜ = view(β, :, t) - Bβₜ₊₁ = view(Bβ, :, t + 1) - predict_previous_state!(βₜ, hmm, Bβₜ₊₁, control_seq[t]) - lmul!(c[t], βₜ) - end - Bβ[:, t1] .= view(B, :, t1) .* view(β, :, t1) - - # State marginals - γ[:, t1:t2] .= view(α, :, t1:t2) .* view(β, :, t1:t2) ./ view(c, t1:t2)' - - # Transition marginals - if transition_marginals - for t in t1:(t2 - 1) - trans = transition_matrix(hmm, control_seq[t]) - mul_rows_cols!(ξ[t], view(α, :, t), trans, view(Bβ, :, t + 1)) - end - ξ[t2] .= zero(R) - end - - return nothing - end - - """ - $(SIGNATURES) - """ - function forward_backward!( - storage::ForwardBackwardStorage, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - transition_marginals::Bool=true, - ) - if seq_ends isa NTuple{1} - for k in eachindex(seq_ends) - _forward_backward!( - storage, hmm, obs_seq, control_seq, seq_ends, k; transition_marginals - ) - end - else - @threads for k in eachindex(seq_ends) - _forward_backward!( - storage, hmm, obs_seq, control_seq, seq_ends, k; transition_marginals - ) - end - end - return nothing - end - - """ - $(SIGNATURES) - - Apply the forward-backward algorithm to infer the posterior state and transition marginals during sequence `obs_seq` for `hmm`. - - Return a tuple `(storage.γ, storage.logL)` where `storage` is of type [`ForwardBackwardStorage`](@ref). - """ - function forward_backward( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector=Fill(nothing, length(obs_seq)); - seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),), - ) - transition_marginals = false - storage = initialize_forward_backward( - hmm, obs_seq, control_seq; seq_ends, transition_marginals - ) - forward_backward!(storage, hmm, obs_seq, control_seq; seq_ends, transition_marginals) - return storage.γ, storage.logL - end - - nothing -end - -# ╔═╡ 7b2b0560-726d-4bc0-acfa-16361da4f0b2 -begin - """ - $(TYPEDEF) - - # Fields - - Only the fields with a description are part of the public API. - - $(TYPEDFIELDS) - """ - struct ViterbiStorage{R} - "most likely state sequence `q[t] = argmaxᵢ ℙ(X[t]=i | Y[1:T])`" - q::Vector{Int} - "one joint loglikelihood per pair of observation sequence and most likely state sequence" - logL::Vector{R} - logB::Matrix{R} - ϕ::Matrix{R} - ψ::Matrix{Int} - end - - """ - $(SIGNATURES) - """ - function initialize_viterbi( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - ) - N, T, K = length(hmm, control_seq[1]), length(obs_seq), length(seq_ends) - R = eltype(hmm, obs_seq[1], control_seq[1]) - q = Vector{Int}(undef, T) - logL = Vector{R}(undef, K) - logB = Matrix{R}(undef, N, T) - ϕ = Matrix{R}(undef, N, T) - ψ = Matrix{Int}(undef, N, T) - return ViterbiStorage(q, logL, logB, ϕ, ψ) - end - - function _viterbi!( - storage::ViterbiStorage{R}, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector, - seq_ends::AbstractVectorOrNTuple{Int}, - k::Integer, - ) where {R} - (; q, logB, ϕ, ψ, logL) = storage - t1, t2 = seq_limits(seq_ends, k) - - logBₜ₁ = view(logB, :, t1) - obs_logdensities!(logBₜ₁, hmm, obs_seq[t1], control_seq[t1], missing) - loginit = log_initialization(hmm, control_seq[t1]) - ϕ[:, t1] .= loginit .+ logBₜ₁ - - for t in (t1 + 1):t2 - logBₜ = view(logB, :, t) - obs_logdensities!( - logBₜ, hmm, obs_seq[t], control_seq[t], previous_obs(hmm, obs_seq, t) - ) - logtrans = log_transition_matrix(hmm, control_seq[t - 1]) - ϕₜ, ϕₜ₋₁ = view(ϕ, :, t), view(ϕ, :, t - 1) - ψₜ = view(ψ, :, t) - argmaxplus_transmul!(ϕₜ, ψₜ, logtrans, ϕₜ₋₁) - ϕₜ .+= logBₜ - end - - ϕₜ₂ = view(ϕ, :, t2) - q[t2] = argmax(ϕₜ₂) - logL[k] = ϕ[q[t2], t2] - for t in (t2 - 1):-1:t1 - q[t] = ψ[q[t + 1], t + 1] - end - - @argcheck isfinite(logL[k]) - return nothing - end - - """ - $(SIGNATURES) - """ - function viterbi!( - storage::ViterbiStorage{R}, - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector; - seq_ends::AbstractVectorOrNTuple{Int}, - ) where {R} - if seq_ends isa NTuple{1} - for k in eachindex(seq_ends) - _viterbi!(storage, hmm, obs_seq, control_seq, seq_ends, k) - end - else - @threads for k in eachindex(seq_ends) - _viterbi!(storage, hmm, obs_seq, control_seq, seq_ends, k) - end - end - return nothing - end - - """ - $(SIGNATURES) - - Apply the Viterbi algorithm to infer the most likely state sequence corresponding to `obs_seq` for `hmm`. - - Return a tuple `(storage.q, storage.logL)` where `storage` is of type [`ViterbiStorage`](@ref). - """ - function viterbi( - hmm::AbstractHMM, - obs_seq::AbstractVector, - control_seq::AbstractVector=Fill(nothing, length(obs_seq)); - seq_ends::AbstractVectorOrNTuple{Int}=(length(obs_seq),), - ) - storage = initialize_viterbi(hmm, obs_seq, control_seq; seq_ends) - viterbi!(storage, hmm, obs_seq, control_seq; seq_ends) - return storage.q, storage.logL - end - - nothing -end - -# ╔═╡ 519901db-5fb9-4323-822a-ff04e28285c9 -forward(hmm₀, true_obs_seq) - -# ╔═╡ 24b2383c-cf58-46a8-84f7-d00ddc1bbbe4 -forward_backward(hmm₀, true_obs_seq) - -# ╔═╡ be40ed94-e6e2-4838-98db-83d1de9eb830 -viterbi(hmm₀, true_obs_seq) - -# ╔═╡ 19b64180-184c-435f-9be6-40dd9f482f78 -pᵏᵢⱼ = [[0.4 0.4 0.2; 0.3 0.3 0.4], [0.25 0.5 0.25; 0.7 0.15 0.15]] - -# ╔═╡ 1bb32d8a-1961-4709-93ff-aa6f988f0ce5 -hmm = DiscreteCARHMM(uᵏᵢ, tᵏᵢⱼ, eᵏʳᵢⱼ, pᵏᵢⱼ) - -# ╔═╡ 3968aa53-4022-4299-b3ea-5d4086872935 -forward(hmm, true_obs_seq, true_control_seq) - -# ╔═╡ 48a686a4-a996-4c1b-8c48-f72d2b8aff4f -forward_backward(hmm, true_obs_seq, true_control_seq) - -# ╔═╡ 7d6012de-53b6-491d-95e4-eff55c2b80b5 -viterbi(hmm, true_obs_seq, true_control_seq) - -# ╔═╡ 3b034316-b7c8-48cc-911a-fdcea177cf2e -md""" -## 9. Comparison -""" - -# ╔═╡ 05e3d2f4-f793-44ce-a9f8-23e02b6dc303 -#https://cran.r-project.org/web/packages/LMest/vignettes/vignetteLMest.html - -# ╔═╡ f3bc9ea7-9e41-47f4-9d3d-fbc7d35bf5a2 -#https://github.com/TheoMichelot/hmmTMB - -# ╔═╡ 78c2591b-063a-4917-936d-28f17686eb7b -#https://github.com/hmmlearn/hmmlearn - -# ╔═╡ 00000000-0000-0000-0000-000000000001 -PLUTO_PROJECT_TOML_CONTENTS = """ -[deps] -ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" -ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" - -[compat] -ArgCheck = "~2.3.0" -ChainRulesCore = "~1.25.0" -DensityInterface = "~0.4.0" -Distributions = "~0.25.113" -DocStringExtensions = "~0.9.3" -FillArrays = "~1.13.0" -PlutoUI = "~0.7.60" -StatsAPI = "~1.7.0" -StatsFuns = "~1.3.2" -""" - -# ╔═╡ 00000000-0000-0000-0000-000000000002 -PLUTO_MANIFEST_TOML_CONTENTS = """ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.11.1" -manifest_format = "2.0" -project_hash = "cb45b4653fe301fab6134ac7f13a025a6c4ad1f3" - -[[deps.AbstractPlutoDingetjes]] -deps = ["Pkg"] -git-tree-sha1 = "6e1d2a35f2f90a4bc7c2ed98079b2ba09c35b83a" -uuid = "6e696c72-6542-2067-7265-42206c756150" -version = "1.3.2" - -[[deps.AliasTables]] -deps = ["PtrArrays", "Random"] -git-tree-sha1 = "9876e1e164b144ca45e9e3198d0b689cadfed9ff" -uuid = "66dad0bd-aa9a-41b7-9441-69ab47430ed8" -version = "1.1.3" - -[[deps.ArgCheck]] -git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4" -uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197" -version = "2.3.0" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.2" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -version = "1.11.0" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" -version = "1.11.0" - -[[deps.ChainRulesCore]] -deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "3e4b134270b372f2ed4d4d0e936aabaefc1802bc" -uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.25.0" -weakdeps = ["SparseArrays"] - - [deps.ChainRulesCore.extensions] - ChainRulesCoreSparseArraysExt = "SparseArrays" - -[[deps.ColorTypes]] -deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d" -uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.11.5" - -[[deps.Compat]] -deps = ["TOML", "UUIDs"] -git-tree-sha1 = "8ae8d32e09f0dcf42a36b90d4e17f5dd2e4c4215" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.16.0" -weakdeps = ["Dates", "LinearAlgebra"] - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - -[[deps.CompilerSupportLibraries_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.1+0" - -[[deps.DataAPI]] -git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" -uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.16.0" - -[[deps.DataStructures]] -deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.20" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" -version = "1.11.0" - -[[deps.DensityInterface]] -deps = ["InverseFunctions", "Test"] -git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b" -uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d" -version = "0.4.0" - -[[deps.Distributions]] -deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "3101c32aab536e7a27b1763c0797dba151b899ad" -uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.113" -weakdeps = ["ChainRulesCore", "DensityInterface", "Test"] - - [deps.Distributions.extensions] - DistributionsChainRulesCoreExt = "ChainRulesCore" - DistributionsDensityInterfaceExt = "DensityInterface" - DistributionsTestExt = "Test" - -[[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" -uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" -version = "1.11.0" - -[[deps.FillArrays]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "6a70198746448456524cb442b8af316927ff3e1a" -uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "1.13.0" -weakdeps = ["PDMats", "SparseArrays", "Statistics"] - - [deps.FillArrays.extensions] - FillArraysPDMatsExt = "PDMats" - FillArraysSparseArraysExt = "SparseArrays" - FillArraysStatisticsExt = "Statistics" - -[[deps.FixedPointNumbers]] -deps = ["Statistics"] -git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172" -uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" -version = "0.8.5" - -[[deps.HypergeometricFunctions]] -deps = ["LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] -git-tree-sha1 = "b1c2585431c382e3fe5805874bda6aea90a95de9" -uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.25" - -[[deps.Hyperscript]] -deps = ["Test"] -git-tree-sha1 = "179267cfa5e712760cd43dcae385d7ea90cc25a4" -uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91" -version = "0.0.5" - -[[deps.HypertextLiteral]] -deps = ["Tricks"] -git-tree-sha1 = "7134810b1afce04bbc1045ca1985fbe81ce17653" -uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2" -version = "0.9.5" - -[[deps.IOCapture]] -deps = ["Logging", "Random"] -git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770" -uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.5" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -version = "1.11.0" - -[[deps.InverseFunctions]] -git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb" -uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.17" -weakdeps = ["Dates", "Test"] - - [deps.InverseFunctions.extensions] - InverseFunctionsDatesExt = "Dates" - InverseFunctionsTestExt = "Test" - -[[deps.IrrationalConstants]] -git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" -uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.2.2" - -[[deps.JLLWrappers]] -deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b" -uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.6.1" - -[[deps.JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" -uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.4" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.6.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -version = "1.11.0" - -[[deps.LibGit2_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] -uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.7.2+0" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -version = "1.11.0" - -[[deps.LinearAlgebra]] -deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] -uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -version = "1.11.0" - -[[deps.LogExpFunctions]] -deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea" -uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.28" - - [deps.LogExpFunctions.extensions] - LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" - LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" - LogExpFunctionsInverseFunctionsExt = "InverseFunctions" - - [deps.LogExpFunctions.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" - InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -version = "1.11.0" - -[[deps.MIMEs]] -git-tree-sha1 = "65f28ad4b594aebe22157d6fac869786a255b7eb" -uuid = "6c6e2e6c-3030-632d-7369-2d6c69616d65" -version = "0.1.4" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -version = "1.11.0" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.6+0" - -[[deps.Missings]] -deps = ["DataAPI"] -git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" -uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.2.0" - -[[deps.Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" -version = "1.11.0" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.12.12" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.OpenBLAS_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] -uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.27+1" - -[[deps.OpenLibm_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+2" - -[[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" -uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" - -[[deps.OrderedCollections]] -git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.3" - -[[deps.PDMats]] -deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65" -uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.31" - -[[deps.Parsers]] -deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" -uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.8.1" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.11.0" - - [deps.Pkg.extensions] - REPLExt = "REPL" - - [deps.Pkg.weakdeps] - REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.PlutoUI]] -deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "FixedPointNumbers", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "MIMEs", "Markdown", "Random", "Reexport", "URIs", "UUIDs"] -git-tree-sha1 = "eba4810d5e6a01f612b948c9fa94f905b49087b0" -uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8" -version = "0.7.60" - -[[deps.PrecompileTools]] -deps = ["Preferences"] -git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" -uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.2.1" - -[[deps.Preferences]] -deps = ["TOML"] -git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6" -uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.4.3" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -version = "1.11.0" - -[[deps.PtrArrays]] -git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f" -uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d" -version = "1.2.1" - -[[deps.QuadGK]] -deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "cda3b045cf9ef07a08ad46731f5a3165e56cf3da" -uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.11.1" - - [deps.QuadGK.extensions] - QuadGKEnzymeExt = "Enzyme" - - [deps.QuadGK.weakdeps] - Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" - -[[deps.Random]] -deps = ["SHA"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -version = "1.11.0" - -[[deps.Reexport]] -git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" -uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "1.2.2" - -[[deps.Rmath]] -deps = ["Random", "Rmath_jll"] -git-tree-sha1 = "852bd0f55565a9e973fcfee83a84413270224dc4" -uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" -version = "0.8.0" - -[[deps.Rmath_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "58cdd8fb2201a6267e1db87ff148dd6c1dbd8ad8" -uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" -version = "0.5.1+0" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -version = "1.11.0" - -[[deps.SortingAlgorithms]] -deps = ["DataStructures"] -git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" -uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.2.1" - -[[deps.SparseArrays]] -deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] -uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.11.0" - -[[deps.SpecialFunctions]] -deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14" -uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.4.0" -weakdeps = ["ChainRulesCore"] - - [deps.SpecialFunctions.extensions] - SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" - -[[deps.Statistics]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0" -uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.11.1" -weakdeps = ["SparseArrays"] - - [deps.Statistics.extensions] - SparseArraysExt = ["SparseArrays"] - -[[deps.StatsAPI]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" -uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" - -[[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "5cf7606d6cef84b543b483848d4ae08ad9832b21" -uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.3" - -[[deps.StatsFuns]] -deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] -git-tree-sha1 = "b423576adc27097764a90e163157bcfc9acf0f46" -uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "1.3.2" -weakdeps = ["ChainRulesCore", "InverseFunctions"] - - [deps.StatsFuns.extensions] - StatsFunsChainRulesCoreExt = "ChainRulesCore" - StatsFunsInverseFunctionsExt = "InverseFunctions" - -[[deps.SuiteSparse]] -deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] -uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" - -[[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] -uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "7.7.0+0" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -version = "1.11.0" - -[[deps.Tricks]] -git-tree-sha1 = "7822b97e99a1672bfb1b49b668a6d46d58d8cbcb" -uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775" -version = "0.1.9" - -[[deps.URIs]] -git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b" -uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.5.1" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -version = "1.11.0" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" -version = "1.11.0" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+1" - -[[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.11.0+0" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.59.0+0" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+2" -""" - -# ╔═╡ Cell order: -# ╠═455866a4-a5ef-11ef-04ff-f347c85da0d8 -# ╠═b8b4c6ec-a8f3-4609-a0c0-15652c0d15be -# ╟─a5c89117-9642-4c28-bf5e-7c93184934b4 -# ╠═6444a4ab-5992-4a39-a2e9-8301f9a6ba91 -# ╟─5d513c86-0595-4fee-8312-9927ea731fde -# ╠═4d76c441-429f-4500-9c26-6247abad17ed -# ╟─b6abd844-345b-4c36-b6c0-699fa5c45fd9 -# ╠═dc54a9a0-4fd6-4f64-9067-775d7674ebc1 -# ╟─ebfa27fc-bdb0-4aaa-8930-85025df33b7c -# ╠═f4fe5934-b4a3-412c-89e2-688bdffd7dd2 -# ╠═c90ce230-7640-4337-b6ae-e9369564e8e7 -# ╟─f2599655-855a-4c41-a6ed-c22c16d792a1 -# ╠═4bcf1079-13c6-418e-a37e-ce71e5cc74f1 -# ╟─e35f3103-c5ec-4a6f-bfc0-76326f8eed54 -# ╠═0a903f68-79dc-4a45-b0fe-e5e17eff3b6f -# ╟─71193717-7f4e-4db6-9669-149e7d6405af -# ╠═7b2b0560-726d-4bc0-acfa-16361da4f0b2 -# ╟─ffb943ac-b403-4f71-a19b-0875748674ed -# ╠═3a00c710-e6e6-4b5e-875e-7bbcb173d8a3 -# ╠═e7acfa37-371b-437c-aa78-47e37500d414 -# ╠═8d99e582-fed1-42bd-8196-5f73cffe5a05 -# ╠═c3ed40cc-0921-411e-aa8b-547d35fa1842 -# ╠═c24abed8-65bc-4bb8-9444-64377c90be21 -# ╠═4b6ecf0d-6330-4fc5-b31a-f49bbc03d397 -# ╠═1c825864-2332-4e99-8d3a-4583f4fb56e2 -# ╠═5e9a1a49-0356-48ec-ad85-ed5c49653d2d -# ╠═519901db-5fb9-4323-822a-ff04e28285c9 -# ╠═24b2383c-cf58-46a8-84f7-d00ddc1bbbe4 -# ╠═be40ed94-e6e2-4838-98db-83d1de9eb830 -# ╠═3b669bd8-21e1-42d5-b2d6-45340c57e1e2 -# ╠═dfa9f848-b928-4299-8347-2ac492958df7 -# ╠═b454b317-29a9-41ec-82a7-9c48a88f2576 -# ╠═b1377b4b-ce15-4eb8-83cb-740a58673ac6 -# ╠═19b64180-184c-435f-9be6-40dd9f482f78 -# ╠═1bb32d8a-1961-4709-93ff-aa6f988f0ce5 -# ╠═3968aa53-4022-4299-b3ea-5d4086872935 -# ╠═48a686a4-a996-4c1b-8c48-f72d2b8aff4f -# ╠═7d6012de-53b6-491d-95e4-eff55c2b80b5 -# ╟─3b034316-b7c8-48cc-911a-fdcea177cf2e -# ╠═05e3d2f4-f793-44ce-a9f8-23e02b6dc303 -# ╠═f3bc9ea7-9e41-47f4-9d3d-fbc7d35bf5a2 -# ╠═78c2591b-063a-4917-936d-28f17686eb7b -# ╟─00000000-0000-0000-0000-000000000001 -# ╟─00000000-0000-0000-0000-000000000002