PF: two-filter smoother

pf_smoother_two_filter() function implements the two-filter particle smoother (Fearnhead et al., 2010).

Usage

set_vmap(.map = NULL, .mobility = NULL, .vmap = NULL, .plot = FALSE, ...)

pf_smoother_two_filter(
  .n_particle = NULL,
  .n_sim = 100L,
  .cache = TRUE,
  .batch = NULL,
  .collect = TRUE,
  .progress = julia_progress(),
  .verbose = getOption("patter.verbose")
)

Arguments

.map, .mobility, .vmap, .plot, ...

(optional) 'Validity map' arguments for set_map(), used for two-dimensional states.

• .map is a terra::SpatRaster that defines the study area of the simulation (see pf_filter()). On Linux, this argument can only be used safely if JULIA_SESSION = "FALSE".

• .mobility is a numeric value that defines the maximum moveable distance between two time steps (e.g., .timeline[1] and .timeline[2] in pf_filter()).

• .vmap is a terra::SpatRaster (supported on Windows or MacOS), or a file path to a raster (supported on MacOS, Windows and Linux), that defines the validity map (see set_map()). This can be supplied, from a previous implementation of set_vmap() or the internal function spatVmap(), instead of .map and .mobility to avoid re-computation.

• .plot is a logical variable that defines whether or not to plot the map.

• ... is a placeholder for additional arguments, passed to terra::plot(), if .plot = TRUE.

The validity map is only set in Julia if JULIA_SESSION = "TRUE".

On Linux, the validity map cannot be created and set in the same R session. Running the function with .map and .mobility will create, but not set, the map. Write the map to file and then rerun the function with .vmap specified to set the map safely in Julia.

.n_particle

(optional) An integer that defines the number of particles to smooth.

• If specified, a sub-sample of .n_particles is used.

• Otherwise, .n_particle = NULL uses all particles from the filter.

.n_sim

An integer that defines the number of Monte Carlo simulations.

.cache

A logical variable that defines whether or not to pre-compute and cache movement-density normalisation constants for each unique particle.

.batch

(optional) Batching controls:

• If pf_filter() was implemented with .batch = NULL, leave .batch = NULL here.

• Otherwise, .batch must be specified. Pass a character vector of .jld2 file paths to write particles in sequential batches to file (as Julia Matrix{<:State} objects) to .batch; for example: ./smo-1.jld2, ./smo-2.jld2, .... You must use the same number of batches as in pf_filter(). .batch is implemented as in that function.

.collect

A logical variable that defines whether or not to collect outputs from the Julia session in R.

.progress

Progress controls (see patter-progress for supported options). If enabled, one progress bar is shown for each .batch.

.verbose

User output control (see patter-progress for supported options).

Value

• set_vmap():

  • set_vmap() returns the validity map (a terra::SpatRaster), invisibly;

• pf_smoother_two_filter():

  • Patter.particle_smoother_two_filter() creates a NamedTuple in the Julia session (named ptf). If .batch = NULL, the NamedTuple contains particles (states) ; otherwise, the states element is nothing and states are written to .jld2 files (as a variable named xsmo). If .collect = TRUE, pf_smoother_two_filter() collects the outputs in R as a pf_particles object (the states element is NULL is .batch is used). Otherwise, invisible(NULL) is returned.

Details

The two-filter smoother smooths particle samples from the particle filter (pf_filter()). Particles from a forward and backward filter run are required in the Julia workspace (as defined by pf_filter()). The backend function Patter.particle_smoother_two_filter() does the work. Essentially, the function runs a simulation backwards in time and re-samples particles in line with the probability density of movements between each combination of states from the backward filter at time t and states from the forward filter at time t - 1. The time complexity of the algorithm is thus \(O(TN^2)\). The probability density of movements is evaluated by Patter.logpdf_step() and Patter.logpdf_move(). If individual states are two-dimensional (see StateXY), a validity map can be pre-defined in Julia via set_vmap() to speed up probability calculations. The validity map is defined as the set of valid (non-NA or bordering) locations on the .map, shrunk by .mobility. Within this region, the probability density of movement between two states can be calculated directly. Otherwise, a Monte Carlo simulation, of .n_sim iterations, is required to compute the normalisation constant (accounting for movements into inhospitable areas, or beyond the boundaries of the study area). For movement models for which the density only depends on the particle states, set .cache = TRUE to pre-compute and cache normalisation constants for improved speed.

References

Fearnhead, P. et al. (2010). A sequential smoothing algorithm with linear computational cost. Biometrika 97, 447–464. https://doi.org/10.1093/biomet/asq013.

See also

Particle filters and smoothers sample states (particles) that represent the possible locations of an individual through time, accounting for all data and the individual's movement.

• To simulate artificial datasets, see sim_*() functions (especially sim_path_walk(), sim_array() and sim_observations()).

• To assemble real-world datasets for the filter, see assemble_*() functions.

• pf_filter() runs the filter:

  • For state types, see State;

  • For observation models, see ModelObs;

  • For movement models, see ModelMove;

• To run particle smoothing, use pf_smoother_two_filter().

• To map emergent patterns of space use, use a map_*() function (such as map_pou(), map_dens() and map_hr()).

Author

Edward Lavender

Examples

if (patter_run()) {

  library(JuliaCall)

  #### Set up example
  # Set up the particle filter with an example dataset
  # (See `?pf_filter()` for the full workflow)
  setup <- example_setup("pf_smoother_two_filter")
  map   <- setup$map
  args  <- setup$pf_filter_args
  # Run the particle filter forwards
  args$.direction <- "forward"
  fwd <- do.call(pf_filter, args)
  # Run the particle filter backwards
  args$.direction <- "backward"
  bwd <- do.call(pf_filter, args)

  #### Example (1): Implement the smoother with default options
  # Run the smoother
  # * This uses objects defined by `pf_filter()` in `Julia`
  # (set_vmap() is explained below)
  smo <- pf_smoother_two_filter()
  # The filter returns a `pf_particles`-class object
  # (See `?pf_filter()` for examples)
  class(smo)
  summary(smo)

  #### Example (2): Implement the smoother using a validity map
  # We can use a validity map b/c `.state` = "StateXY"
  args$.state
  # To define the validity map, define:
  # * `.map`
  # * `.mobility`, which we can see here is 750 m:
  args$.model_move
  # Run the smoother
  set_vmap(.map = map, .mobility = 750.0, .plot = TRUE)
  smo <- pf_smoother_two_filter()
  # Reset `vmap` in `Julia` to run the smoother for other state types
  # ... in the same R session:
  set_vmap()

  #### Example (3): Implement the smoother with a sub-sample of particles
  # This is useful for quick tests
  set_vmap(.map = map, .mobility = 750.0)
  smo <- pf_smoother_two_filter(.n_particle = 50L)

  #### Example (4): Adjust the number of MC simulations
  # set_vmap(.map = map, .mobility = 750.0)
  smo <- pf_smoother_two_filter(.n_sim = 1000L)

  #### Example (5): Batching workflow for filtering and smoothing
  # (Use `.batch` to reduce memory demand)

  folder <- file.path(tempdir(), "batch-outputs")
  dir.create(folder)

  ## Run forward filter with batching
  batch_fwd       <- file.path(folder, c("fwd-1.jld2", "fwd-2.jld2", "fwd-3.jld2"))
  args$.direction <- "forward"
  args$.batch     <- batch_fwd
  args$.collect   <- TRUE
  fwd <- do.call(pf_filter, args)
  # In the output, the 'states' element is null
  fwd$states
  # Other elements are as expected
  fwd$callstats
  head(fwd$diagnostics)
  # Confirm that batch files exist
  stopifnot(all(file.exists(batch_fwd)))

  ## Run backward filter with batching:
  batch_bwd       <- file.path(folder, c("bwd-1.jld2", "bwd-2.jld2", "bwd-3.jld2"))
  args$.direction <- "backward"
  args$.batch     <- batch_bwd
  bwd <- do.call(pf_filter, args)
  summary(bwd)
  stopifnot(all(file.exists(batch_bwd)))

  ## Run smoothing with batching
  # set_vmap(.map = map, .mobility = 750.0)
  batch_smo <- file.path(folder, c("smo-1.jld2", "smo-2.jld2", "smo-3.jld2"))
  smo       <- pf_smoother_two_filter(.batch = batch_smo)
  summary(smo)
  stopifnot(all(file.exists(batch_smo)))

  ## Collate outputs in R
  julia_command('
    smo_states = hcat([f["xsmo"] for f in map(jldopen, batch_smo)]...);')
  smo$states <- julia_eval('
    Patter.r_get_states(smo_states, collect(1:length(timeline)), timeline);
    ')
  head(smo$states)

  #### Example (6): Analyse smoothed particles
  # * See `map_*()` functions (e.g., `?map_dens()`) to map utilisation distributions


  # Cleanup
  file_cleanup(folder)
  set_vmap()
}
#> `patter::julia_connect()` called @ 2025-04-22 09:31:41... 
#> ... Running `Julia` setup via `JuliaCall::julia_setup()`... 
#> ... Validating Julia installation... 
#> ... Setting up Julia project... 
#> ... Handling dependencies... 
#> ... `Julia` set up with 11 thread(s). 
#> `patter::julia_connect()` call ended @ 2025-04-22 09:31:41 (duration: ~0 sec(s)). 


#> `patter::pf_filter()` called @ 2025-04-22 09:31:41... 
#> `patter::pf_filter_init()` called @ 2025-04-22 09:31:41... 
#> ... 09:31:41: Setting initial states... 
#> ... 09:31:41: Setting observations dictionary... 
#> `patter::pf_filter_init()` call ended @ 2025-04-22 09:31:41 (duration: ~0 sec(s)). 
#> ... 09:31:41: Running filter... 
#> Message: On iteration 1 ...
#>   
#> Message: Running filter for batch 1 / 1 ...
#>   
#> ... 09:31:42: Collating outputs... 
#> `patter::pf_filter()` call ended @ 2025-04-22 09:31:42 (duration: ~1 sec(s)). 
#> `patter::pf_filter()` called @ 2025-04-22 09:31:42... 
#> `patter::pf_filter_init()` called @ 2025-04-22 09:31:42... 
#> ... 09:31:42: Setting initial states... 
#> ... 09:31:42: Setting observations dictionary... 
#> `patter::pf_filter_init()` call ended @ 2025-04-22 09:31:42 (duration: ~0 sec(s)). 
#> ... 09:31:42: Running filter... 
#> Message: On iteration 1 ...
#>   
#> Message: Running filter for batch 1 / 1 ...
#>   
#> ... 09:31:42: Collating outputs... 
#> `patter::pf_filter()` call ended @ 2025-04-22 09:31:42 (duration: ~0 sec(s)). 
#> `patter::pf_smoother_two_filter()` called @ 2025-04-22 09:31:42... 
#> ... 09:31:42: Running smoother... 
#> Message: Running smoother for batch 1 / 1 ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> ... 09:31:43: Collating outputs... 
#> `patter::pf_smoother_two_filter()` call ended @ 2025-04-22 09:31:43 (duration: ~1 sec(s)). 

#> `patter::pf_smoother_two_filter()` called @ 2025-04-22 09:31:43... 
#> ... 09:31:43: Running smoother... 
#> Message: Running smoother for batch 1 / 1 ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> ... 09:31:44: Collating outputs... 
#> `patter::pf_smoother_two_filter()` call ended @ 2025-04-22 09:31:44 (duration: ~1 sec(s)). 
#> `patter::pf_smoother_two_filter()` called @ 2025-04-22 09:31:44... 
#> ... 09:31:44: Running smoother... 
#> Message: Running smoother for batch 1 / 1 ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> ... 09:31:44: Collating outputs... 
#> `patter::pf_smoother_two_filter()` call ended @ 2025-04-22 09:31:44 (duration: ~0 sec(s)). 
#> `patter::pf_smoother_two_filter()` called @ 2025-04-22 09:31:44... 
#> ... 09:31:44: Running smoother... 
#> Message: Running smoother for batch 1 / 1 ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> ... 09:31:44: Collating outputs... 
#> `patter::pf_smoother_two_filter()` call ended @ 2025-04-22 09:31:44 (duration: ~0 sec(s)). 
#> `patter::pf_filter()` called @ 2025-04-22 09:31:44... 
#> `patter::pf_filter_init()` called @ 2025-04-22 09:31:44... 
#> ... 09:31:44: Setting initial states... 
#> ... 09:31:44: Setting observations dictionary... 
#> `patter::pf_filter_init()` call ended @ 2025-04-22 09:31:44 (duration: ~0 sec(s)). 
#> ... 09:31:44: Running filter... 
#> Message: On iteration 1 ...
#>   
#> Message: Running filter for batch 1 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> Message: Running filter for batch 2 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> Message: Running filter for batch 3 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> ... 09:31:46: Collating outputs... 
#> `patter::pf_filter()` call ended @ 2025-04-22 09:31:46 (duration: ~2 sec(s)). 
#> `patter::pf_filter()` called @ 2025-04-22 09:31:46... 
#> `patter::pf_filter_init()` called @ 2025-04-22 09:31:46... 
#> ... 09:31:46: Setting initial states... 
#> ... 09:31:46: Setting observations dictionary... 
#> `patter::pf_filter_init()` call ended @ 2025-04-22 09:31:46 (duration: ~0 sec(s)). 
#> ... 09:31:46: Running filter... 
#> Message: On iteration 1 ...
#>   
#> Message: Running filter for batch 1 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> Message: Running filter for batch 2 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> Message: Running filter for batch 3 / 3 ...
#>   
#> Message: Writing particles to disk ...
#>   
#> ... 09:31:46: Collating outputs... 
#> `patter::pf_filter()` call ended @ 2025-04-22 09:31:46 (duration: ~0 sec(s)). 
#> `patter::pf_smoother_two_filter()` called @ 2025-04-22 09:31:46... 
#> ... 09:31:46: Running smoother... 
#> Message: Running smoother for batch 1 / 3 ...
#>   
#> Message: Reading particles from disk ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> Message: Writing smoothed particles to disk ...
#>   
#> Message: Running smoother for batch 2 / 3 ...
#>   
#> Message: Reading particles from disk ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> Message: Writing smoothed particles to disk ...
#>   
#> Message: Running smoother for batch 3 / 3 ...
#>   
#> Message: Reading particles from disk ...
#>   
#> Message: Precomputing normalisation constants for batch ...
#>   
#> Message: Smoothing ...
#>   
#> Message: Writing smoothed particles to disk ...
#>   
#> ... 09:31:47: Collating outputs... 
#> `patter::pf_smoother_two_filter()` call ended @ 2025-04-22 09:31:47 (duration: ~1 sec(s)).