Skip to content
Snippets Groups Projects
Commit a7c6db10 authored by Rob Moss's avatar Rob Moss
Browse files

Avoid invalid assumptions about particle ordering

When starting from a cached state, we cannot assume that the initial
particle ordering is consistent with previous time-steps.
parent da624c09
No related branches found
No related tags found
No related merge requests found
......@@ -30,7 +30,7 @@ def expect(ctx, time, unit, period, prev, curr):
raise ValueError("Unknown observation type '{}'".format(unit))
def log_llhd_of(ctx, hist, hist_ix, obs, max_back=None):
def log_llhd_of(ctx, hist, hist_ix, obs, max_back=0):
"""Return the log-likelihood of obtaining observations from each particle.
:param ctx: The simulation context.
......@@ -39,7 +39,7 @@ def log_llhd_of(ctx, hist, hist_ix, obs, max_back=None):
:param obs: The observation(s) that have been made.
:param max_back: The number of time-steps into the past when the most
recent resampling occurred (i.e., how far back the current particle
ordering is guaranteed to persist; default is ``None``, no limit).
ordering is guaranteed to persist; default is 0, no guarantee).
:returns: An array containing the log-likelihood for each particle.
"""
......@@ -58,7 +58,7 @@ def log_llhd_of(ctx, hist, hist_ix, obs, max_back=None):
def hist_for(period):
"""Return past state vectors in the appropriate order."""
steps_back = steps_per_unit * period
same_ixs = max_back is None or max_back >= steps_back
same_ixs = max_back >= steps_back
if same_ixs:
if steps_back > hist_ix:
# If the observation period starts before the beginning of the
......
......@@ -8,7 +8,7 @@ from . import obs as obs_mod
from . import state as state_mod
def reweight(ctx, hist, hist_ix, obs, max_back=None):
def reweight(ctx, hist, hist_ix, obs, max_back):
"""Adjust particle weights in response to some observation(s).
:param params: The simulation parameters.
......@@ -17,7 +17,7 @@ def reweight(ctx, hist, hist_ix, obs, max_back=None):
:param obs: The observation(s) that have been made.
:param max_back: The number of time-steps into the past when the most
recent resampling occurred (i.e., how far back the current particle
ordering is guaranteed to persist; default is ``None``, no limit).
ordering is guaranteed to persist).
:returns: A tuple; the first element (*bool*) indicates whether resampling
is required, the second element (*float*) is the **effective** number
......@@ -200,7 +200,12 @@ def run(ctx, start, end, streams, state=None,
# The time of the previous time-step (if any).
most_recent = None
# The time-step number of the most recent resampling (if any).
last_rs = None
# NOTE: the first time-step is number 1 and updates hist[1] given hist[0].
# So we set this value to zero to indicate that we can step back as far as
# the beginning of this simulation; this is important when resuming from
# cached states where we cannot assume anything about the history prior to
# the first time-step.
last_rs = 0
# The index of the current time-step in the state matrix.
hist_ix = None
......@@ -237,10 +242,7 @@ def run(ctx, start, end, streams, state=None,
hist[hist_ix, :, :-2] = 0
# Determine how many time-steps back the most recent resampling was.
if last_rs is None:
max_back = None
else:
max_back = (step_num - last_rs)
max_back = (step_num - last_rs)
# Simulate the current time-step.
resampled = step(ctx, hist, hist_ix, step_num, when, obs,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment