From 72b0ecbd0e9c531a2bb488155dc73cc4d1efdb71 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 3 Jun 2024 08:18:10 +0200 Subject: [PATCH 1/2] remove the unused argument `df` --- docs/src/beetle_example.md | 2 +- docs/src/index.md | 2 +- src/PFtypes.jl | 6 ++++-- test/runtests.jl | 4 ++-- 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/docs/src/beetle_example.md b/docs/src/beetle_example.md index 9fa5708f..d3d9ba42 100644 --- a/docs/src/beetle_example.md +++ b/docs/src/beetle_example.md @@ -82,7 +82,7 @@ nothing # hide In this example, we have no control inputs, we thus define a vector of only zeros. We then solve the forward filtering problem and plot the results. ```@example beetle u = zeros(length(y)) -pf = AuxiliaryParticleFilter(AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0)) +pf = AuxiliaryParticleFilter(AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0)) T = length(y) sol = forward_trajectory(pf,u[1:T],y[1:T]) (; x,w,we,ll) = sol diff --git a/docs/src/index.md b/docs/src/index.md index 59d8624d..ca1840e6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -289,7 +289,7 @@ nothing # hide We now create the [`AdvancedParticleFilter`](@ref) and use it in the same way as the other filters: ```@example lingauss -apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0) +apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0) sol = forward_trajectory(apf, u, y, p) ``` diff --git a/src/PFtypes.jl b/src/PFtypes.jl index 1d8ddf5b..6630d00b 100644 --- a/src/PFtypes.jl +++ b/src/PFtypes.jl @@ -175,7 +175,7 @@ end """ - AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, dynamics_density, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...) + AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...) This type represents a standard particle filter but affords extra flexibility compared to the [`ParticleFilter`](@ref) type, e.g., non-additive noise in the dynamics and measurement functions. @@ -186,7 +186,6 @@ See the docs for more information: https://baggepinnen.github.io/LowLevelParticl - `dynamics`: A discrete-time dynamics function `(x, u, p, t, noise=false) -> x⁺`. It's important that the `noise` argument defaults to `false`. - `measurement`: A measurement function `(x, u, p, t, noise=false) -> y`. It's important that the `noise` argument defaults to `false`. - `measurement_likelihood`: A function `(x, u, y, p, t)->logl` to evaluate the log-likelihood of a measurement. -- `dynamics_density`: This field is not used by the advanced filter and can be set to `nothing`. - `initial_density`: The distribution of the initial state. - `threads`: use threads to propagate particles in parallel. Only activate this if your dynamics is thread-safe. `SeeToDee.SimpleColloc` is not thread-safe by default due to the use of internal caches, but `SeeToDee.Rk4` is. """ @@ -202,6 +201,9 @@ function AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Fun initial_density, kwargs...) end +# Constructor without dynamics_density which is not used +AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; kwargs...) = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, Normal(), initial_density; kwargs...) + Base.@propagate_inbounds function measurement_equation!(pf::AbstractParticleFilter, u, y, p, t, w = weights(pf)) g = measurement_likelihood(pf) diff --git a/test/runtests.jl b/test/runtests.jl index 656a0d51..395e4e4f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -292,7 +292,7 @@ end T = 200 # Number of time steps N = 500 Random.seed!(0) - apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0, threads=false) + apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0, threads=false) x,u,y = LowLevelParticleFilters.simulate(apf,T,du) # Simuate trajectory using the model in the filter x,u,y = tosvec.((x,u,y)) @time resapf,ll = mean_trajectory(apf, u, y) @@ -302,7 +302,7 @@ end # With threads Random.seed!(0) - apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0, threads=true) + apf = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0, threads=true) x,u,y = LowLevelParticleFilters.simulate(apf,T,du) # Simuate trajectory using the model in the filter x,u,y = tosvec.((x,u,y)) @time resapf,ll = mean_trajectory(apf, u, y) From e7d10d5ba5106fe47ee5e3125b8b5c92d5587243 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 3 Jun 2024 08:40:25 +0200 Subject: [PATCH 2/2] warn if missing distribution for dynamics density --- docs/src/beetle_example.md | 6 +++--- src/PFtypes.jl | 6 +++--- src/smoothing.jl | 1 + src/utils.jl | 2 ++ 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/docs/src/beetle_example.md b/docs/src/beetle_example.md index d3d9ba42..d29bdcce 100644 --- a/docs/src/beetle_example.md +++ b/docs/src/beetle_example.md @@ -82,7 +82,7 @@ nothing # hide In this example, we have no control inputs, we thus define a vector of only zeros. We then solve the forward filtering problem and plot the results. ```@example beetle u = zeros(length(y)) -pf = AuxiliaryParticleFilter(AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, d0)) +pf = AuxiliaryParticleFilter(AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, df, d0)) T = length(y) sol = forward_trajectory(pf,u[1:T],y[1:T]) (; x,w,we,ll) = sol @@ -108,8 +108,8 @@ plot!(xyt[:,1],xyt[:,2], c=:red, lab="measurement") ``` as well as the angle state variable (we subsample the particles to not get sluggish plots) ```@example beetle -fig2 = scatter(to1series(ϕ.(x)'[:,1:5:end])..., m=(:black, 0.03, 2), lab="", size=(500,300), format=:png) -plot!(identity.(xh[:,4]), lab="Filtered angle", legend=:topleft, ylims=(-30, 70)) +fig2 = scatter(to1series(ϕ.(x)'[:,1:5:end])..., m=(:black, 0.03, 2), lab="", size=(500,300), format=:png, dpi=100) +plot!(identity.(xh[:,4]), lab="Filtered angle", legend=:topleft, ylims=(-30, 70), format=:png) ``` The particle plot above indicate that the posterior is multimodal. This phenomenon arises due to the simple model that uses an angle that is allowed to leave the interval ``0-2\pi` rad. In this example, we are not interested in the angle, but rather when the beetle switches mode. The filtering distribution above gives a hint at when this happens, but we will not plot the mode trajectory until we have explored smoothing as well. diff --git a/src/PFtypes.jl b/src/PFtypes.jl index 6630d00b..21559d27 100644 --- a/src/PFtypes.jl +++ b/src/PFtypes.jl @@ -164,7 +164,7 @@ end dynamics::FT measurement::GT measurement_likelihood::GLT - dynamics_density::FDT = Normal() + dynamics_density::FDT = MissingDistribution() initial_density::IDT resample_threshold::Float64 = 0.5 resampling_strategy::RST = ResampleSystematic @@ -201,8 +201,8 @@ function AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Fun initial_density, kwargs...) end -# Constructor without dynamics_density which is not used -AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; kwargs...) = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, Normal(), initial_density; kwargs...) +# Constructor without dynamics_density which is not used except for in smoothing +AdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, initial_density; kwargs...) = AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, MissingDistribution(), initial_density; kwargs...) Base.@propagate_inbounds function measurement_equation!(pf::AbstractParticleFilter, u, y, p, t, w = weights(pf)) diff --git a/src/smoothing.jl b/src/smoothing.jl index 9bae2d8d..71a301c6 100644 --- a/src/smoothing.jl +++ b/src/smoothing.jl @@ -43,6 +43,7 @@ function smooth(pf::AbstractParticleFilter, xf, wf, wef, ll, M, u, y, p=paramete N = num_particles(pf) f = dynamics(pf) df = dynamics_density(pf) + df isa MissingDistribution && error("In order to perform smoothing, you must provide a noise density for the dynamics. Use the constructor `AdvancedParticleFilter(N, dynamics, measurement, measurement_likelihood, dynamics_density, initial_density; kwargs...)` to do this") @assert M <= N "Must extend cache size of bins and j to allow this" xb = Array{particletype(pf)}(undef,M,T) j = resample(ResampleSystematic, wef[:,T], M) diff --git a/src/utils.jl b/src/utils.jl index a81c9cea..2b7ab2ed 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -115,6 +115,8 @@ Mixed value support indicates that the distribution is a mix of continuous and d """ struct Mixed <: ValueSupport end +struct MissingDistribution <: MultivariateDistribution{Mixed} end + """ TupleProduct(v::NTuple{N,UnivariateDistribution})