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

Minor refactor of pypfilt.resample

parent bd740320
Branches
Tags
No related merge requests found
......@@ -39,9 +39,8 @@ def post_regularise(ctx, px, new_px):
# List parameters in the order they appear in the state vector.
smooth_fields = [n for n in field_names if n in smooth_fields]
num_params = len(smooth_fields)
p_ixs = np.array(range(num_params))
if len(p_ixs) == 0:
if num_params == 0:
logger.debug("Post-RPF: no parameters to resample")
return
......@@ -51,21 +50,20 @@ def post_regularise(ctx, px, new_px):
# Check for parameters that are constant (or nearly so) for all particles.
# These parameters must be ignored or the covariance matrix will not be
# positive definite, and the Cholesky decomposition will fail.
p_range = np.ptp(x[:, p_ixs], axis=0)
p_range = np.ptp(x, axis=0)
toln = ctx.params['resample']['reg_toln']
good = p_range >= toln
if not np.all(good):
bad = np.logical_not(good)
msg = "Post-RPF found {} constant parameter(s) at {}".format(
sum(bad), p_ixs[np.where(bad)])
sum(bad), np.array(smooth_fields)[bad])
logger.debug(msg)
# Update the variables related to these parameters.
smooth_fields = [name for name in np.array(smooth_fields)[good]]
num_params = len(smooth_fields)
p_ixs = np.array(range(num_params))
if len(p_ixs) == 0:
if num_params == 0:
logger.debug("Post-RPF: no non-constant parameters to resample")
return
......@@ -75,8 +73,7 @@ def post_regularise(ctx, px, new_px):
# Use a bandwidth that is half that of the optimal bandwidth for a
# Gaussian kernel (when the underlying density is Gaussian with unit
# covariance), to handle multi-model densities.
npar = len(p_ixs)
h = 0.5 * (4 / (count * (npar + 2))) ** (1 / (npar + 4))
h = 0.5 * (4 / (count * (num_params + 2))) ** (1 / (num_params + 4))
# Calculate the Cholesky decomposition of the parameter covariance
# matrix V, which is used to transform independent normal samples
......@@ -107,7 +104,7 @@ def post_regularise(ctx, px, new_px):
return
# Sample the multivariate normal with covariance V and mean of zero.
std_samples = rnd.normal(size=(npar, count))
std_samples = rnd.normal(size=(num_params, count))
scaled_samples = np.transpose(np.dot(a_mat, h * std_samples))
# Add the sampled noise and clip to respect parameter bounds.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment