Skip to content

Commit

Permalink
Improve performance tips part of the section
Browse files Browse the repository at this point in the history
  • Loading branch information
baggepinnen committed Nov 6, 2024
1 parent 9e0900d commit e6ebaf4
Show file tree
Hide file tree
Showing 4 changed files with 17 additions and 11 deletions.
2 changes: 1 addition & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ makedocs(
"Discretization" => "discretization.md",
"Parameter estimation" => "parameter_estimation.md",
"Benchmark" => "benchmark.md",
"High-performance distributions" => "distributions.md",
"Performance tips" => "distributions.md",
"Tutorials" => [
"Kalman-filter tutorial with LowLevelParticleFilters" => "adaptive_kalmanfilter.md",
"Noise tuning and disturbance modeling for Kalman filtering" => "noisetuning.md",
Expand Down
10 changes: 10 additions & 0 deletions docs/src/distributions.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,13 @@
# Performance tips
Use of [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) is recommended for optimal performance when the state dimension is small, e.g., less than about 10-15 for Kalman filters and less than about 100 for particle filters. In the section [Parameter optimization](https://baggepinnen.github.io/LowLevelParticleFilters.jl/dev/parameter_estimation/#Using-an-optimizer) we demonstrate one workflow that makes use of StaticArrays everywhere it is needed for an [`UnscentedKalmanFilter`](@ref) in order to get a completely allocation free filter. The following arrays must be static for this to hold

- The initial state distribution (the vector and matrix passed to `d0 = MvNormal(μ, Σ)` for Kalman filters). If you are performing parameter optimization with gradients derived using ForwardDiff.jl, these must further have the correct element type. How to achieve this is demonstrated in the liked example above.
- Inputs `u` measured outputs `y`.
- In case of Kalman filters, the dynamic model matrices `A`, `B`, `C`, `D` and the covariance matrices `R1`, `R2`.
- The dynamics functions for [`UnscentedKalmanFilter`](@ref) and particle filters must further return static arrays when passed static arrays as inputs.

## Analysis using JET
All flavors of Kalman filters are analyzed for potential runtime dispatch using [JET.jl](https://github.com/aviatesk/JET.jl). This analysis is performed [in the tests](https://github.com/baggepinnen/LowLevelParticleFilters.jl/blob/master/test/test_jet.jl) and generally requires a completely static filter using static arrays internally. See the tests for an example of how to set a filter up this way.

# High performance Distributions
When `using LowLevelParticleFilters`, a number of methods related to distributions are defined for static arrays, making `logpdf` etc. faster. We also provide a new kind of distribution: `TupleProduct <: MultivariateDistribution` that behaves similarly to the `Product` distribution. The `TupleProduct` however stores the individual distributions in a tuple, has compile-time known length and supports `Mixed <: ValueSupport`, meaning that it can be a product of both `Continuous` and `Discrete` dimensions, something not supported by the standard `Product`. Example
Expand Down
10 changes: 5 additions & 5 deletions docs/src/parameter_estimation.md
Original file line number Diff line number Diff line change
Expand Up @@ -271,9 +271,9 @@ nothing # hide
We then define a nonlinear state estimator, we will use the [`UnscentedKalmanFilter`](@ref), and solve the filtering problem. We start by an initial state estimate ``x_0`` that is slightly off for the parameter ``a_1``
```@example paramest
nx = 5
R1 = SMatrix{nx,nx}(Diagonal([0.1, 0.1, 0.1, 0.1, 0.0001]))
R1 = SMatrix{nx,nx}(Diagonal([0.1, 0.1, 0.1, 0.1, 0.0001])) # Use of StaticArrays is generally good for performance
R2 = SMatrix{ny,ny}(Diagonal((1e-2)^2 * ones(ny)))
x0 = SA[2, 2, 3, 3, 0.02]
x0 = SA[2, 2, 3, 3, 0.02] # The SA prefix makes the array static, which is good for performance
kf = UnscentedKalmanFilter(discrete_dynamics_params, measurement, R1, R2, MvNormal(x0, R1); ny, nu)
Expand Down Expand Up @@ -323,7 +323,7 @@ t = 0:Ts:1000
u1 = vcat.(0.25 .* sign.(sin.(2pi/Tperiod .* (t ./ 40).^2)) .+ 0.25)
u2 = vcat.(0.25 .* sign.(sin.(2pi/Tperiod .* (t ./ 40).^2 .+ pi/2)) .+ 0.25)
u = SVector{nu}.(vcat.(u1,u2))
x0 = SA[2.0,2,3,3]
x0 = SA[2.0,2,3,3] # Initial condition, static array for performance
x = LowLevelParticleFilters.rollout(discrete_dynamics, x0, u, p_true)[1:end-1]
y = measurement.(x, u, 0, 0)
y = [y .+ 0.01 .* randn.() for y in y]
Expand All @@ -335,10 +335,10 @@ plot(
```


This time, we define a cost function for the optimizer to optimize, we'll use the sum of squared errors (`sse`). It's important to define the UKF with an initial state distribution with the same element type as the parameter vector so that automatic differentiation through the state estimator works, hence the explicit casting `T.(x0)` and `T.(R1)`.
This time, we define a cost function for the optimizer to optimize, we'll use the sum of squared errors (`sse`). It's important to define the UKF with an initial state distribution with the same element type as the parameter vector so that automatic differentiation through the state estimator works, hence the explicit casting `T.(x0)` and `T.(R1)`. We also make sure to use StaticArrays for the covariance matrices and the initial condition for performance reasons (optional).
```@example paramest
nx = 4
R1 = SMatrix{nx,nx}(Diagonal([0.1, 0.1, 0.1, 0.1]))
R1 = SMatrix{nx,nx}(Diagonal([0.1, 0.1, 0.1, 0.1])) # Use of StaticArrays is generally good for performance
R2 = SMatrix{ny,ny}(Diagonal((1e-2)^2 * ones(ny)))
x0 = SA[2.0, 2, 3, 3]
Expand Down
6 changes: 1 addition & 5 deletions test/test_jet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,7 @@ skf = SqKalmanFilter(_A, _B, _C, 0, eye(nx), eye(ny), d0)
## Test allocations ============================================================
forward_trajectory(kf, u, y)
a = @allocations forward_trajectory(kf, u, y)
@test a <= 15

@test a <= 15 # Allocations occur when the arrays are allocated for saving the data, the important thing is that the number of allocations do not grow with the length of the trajectory (T = 200)

forward_trajectory(ukf, u, y)
a = @allocations forward_trajectory(ukf, u, y)
Expand All @@ -68,6 +67,3 @@ forward_trajectory(skf, u, y)
a = @allocations forward_trajectory(skf, u, y)

@test a <= 50 # was 7 on julia v1.10.6


## Test differentiability ======================================================

0 comments on commit e6ebaf4

Please sign in to comment.