From eb5fb9f0fcbd5105fbe0b78ff097f62c3e54ef57 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:24:32 -0700 Subject: [PATCH 01/24] docs: correct input indexing in SSM and HMM transition equations --- dynamax/hidden_markov_model/inference.py | 14 +++++++------- dynamax/hidden_markov_model/models/abstractions.py | 2 +- dynamax/hidden_markov_model/parallel_inference.py | 6 +++--- dynamax/ssm.py | 4 ++-- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/dynamax/hidden_markov_model/inference.py b/dynamax/hidden_markov_model/inference.py index 964499a84..5a8a95db2 100644 --- a/dynamax/hidden_markov_model/inference.py +++ b/dynamax/hidden_markov_model/inference.py @@ -127,7 +127,7 @@ def hmm_filter( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -175,7 +175,7 @@ def hmm_backward_filter( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -224,7 +224,7 @@ def hmm_two_filter_smoother( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -280,7 +280,7 @@ def hmm_smoother( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -362,7 +362,7 @@ def hmm_fixed_lag_smoother( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. window_size: the number of future steps to use, $L$ transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -468,7 +468,7 @@ def hmm_posterior_mode( Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. @@ -519,7 +519,7 @@ def hmm_posterior_sample( Args: rng: random number generator initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. transition_fn: function that takes in an integer time index and returns a $K \times K$ transition matrix. diff --git a/dynamax/hidden_markov_model/models/abstractions.py b/dynamax/hidden_markov_model/models/abstractions.py index 1dd6b946b..533b0ffcf 100644 --- a/dynamax/hidden_markov_model/models/abstractions.py +++ b/dynamax/hidden_markov_model/models/abstractions.py @@ -230,7 +230,7 @@ def distribution(self, Args: params: transition parameters state: current latent state - inputs: current inputs + inputs: next inputs Returns: conditional distribution of next state. diff --git a/dynamax/hidden_markov_model/parallel_inference.py b/dynamax/hidden_markov_model/parallel_inference.py index d890f8f35..384c7e299 100644 --- a/dynamax/hidden_markov_model/parallel_inference.py +++ b/dynamax/hidden_markov_model/parallel_inference.py @@ -49,7 +49,7 @@ def hmm_filter(initial_probs: Float[Array, " num_states"], Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. Returns: @@ -114,7 +114,7 @@ def hmm_smoother(initial_probs: Float[Array, " num_states"], Args: initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. Returns: @@ -185,7 +185,7 @@ def hmm_posterior_sample(key: Array, Args: key: random number generator initial_distribution: $p(z_1 \mid u_1, \theta)$ - transition_matrix: $p(z_{t+1} \mid z_t, u_t, \theta)$ + transition_matrix: $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$ log_likelihoods: $p(y_t \mid z_t, u_t, \theta)$ for $t=1,\ldots, T$. Returns: diff --git a/dynamax/ssm.py b/dynamax/ssm.py index 0c778ef85..45f067b97 100644 --- a/dynamax/ssm.py +++ b/dynamax/ssm.py @@ -118,10 +118,10 @@ def transition_distribution( Args: params: model parameters $\theta$ state: current latent state $z_t$ - inputs: current inputs $u_t$ + inputs: next inputs $u_{t+1}$ Returns: - conditional distribution of next latent state $p(z_{t+1} \mid z_t, u_t, \theta)$. + conditional distribution of next latent state $p(z_{t+1} \mid z_t, u_{t+1}, \theta)$. """ raise NotImplementedError From f33b4f715be01ac9245c2ddb90f2cfaf474bee35 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:25:50 -0700 Subject: [PATCH 02/24] test: check serial LGSSM input and parameter indexing, covariance lengths, and forecasts --- .../inference_test_utils.py | 56 +++++++++++++++++++ .../linear_gaussian_ssm/input_timing_test.py | 41 ++++++++++++++ dynamax/linear_gaussian_ssm/models_test.py | 32 ++++++++++- 3 files changed, 128 insertions(+), 1 deletion(-) create mode 100644 dynamax/linear_gaussian_ssm/inference_test_utils.py create mode 100644 dynamax/linear_gaussian_ssm/input_timing_test.py diff --git a/dynamax/linear_gaussian_ssm/inference_test_utils.py b/dynamax/linear_gaussian_ssm/inference_test_utils.py new file mode 100644 index 000000000..fb57bbe02 --- /dev/null +++ b/dynamax/linear_gaussian_ssm/inference_test_utils.py @@ -0,0 +1,56 @@ +"""Exact LGSSM results for input and parameter indexing tests.""" + +import numpy as np +import jax.numpy as jnp + +from dynamax.linear_gaussian_ssm.inference import make_lgssm_params + + +def make_timing_case(num_timesteps=3): + """Construct a one- or three-timestep scalar model and its exact reference results.""" + assert num_timesteps in (1, 3) + inputs = jnp.array([[8.0], [2.0], [4.0]])[:num_timesteps] + scale = inputs / 2 + params = make_lgssm_params( + initial_mean=jnp.zeros(1), initial_cov=jnp.eye(1), + dynamics_weights=scale[:, :, None], + dynamics_input_weights=scale[:, :, None], + dynamics_bias=scale, + # The initial state has no transition, so Q[0] must be unused. + dynamics_cov=scale[:, :, None].at[0].set(jnp.nan), + emissions_weights=jnp.eye(1), emissions_cov=jnp.eye(1), + emissions_input_weights=jnp.array([[0.5]]), + ) + emissions = jnp.array([[4.0], [4.0], [5.0]])[:num_timesteps] + + # Innovations are [0, 0, -13], with variances [2, 5/2, 27/5]. + # The backward smoothing gains are [1/3, 3/11]. + if num_timesteps == 3: + means = np.array([-26/27, 1/9, 146/27]) + covariance = np.array([[10/27, 1/9, 2/27], [1/9, 1/3, 2/9], [2/27, 2/9, 22/27]]) + loglik = -0.5 * (3 * np.log(2 * np.pi) + np.log(27) + 169 / (27/5)) + else: + means, covariance = np.array([0.0]), np.array([[0.5]]) + loglik = -0.5 * np.log(4 * np.pi) + expected = dict( + filtered_means=np.array([0, 3, 146/27])[:num_timesteps, None], + filtered_covariances=np.array([1/2, 3/5, 22/27])[:num_timesteps, None, None], + predicted_means=np.array([3, 16])[:num_timesteps - 1, None], + predicted_covariances=np.array([3/2, 22/5])[:num_timesteps - 1, None, None], + smoothed_means=means[:, None], + smoothed_covariances=covariance.diagonal()[:, None, None], + smoothed_cross_covariances=(covariance.diagonal(1) + means[:-1] * means[1:])[:, None, None], + marginal_loglik=loglik, + joint_covariance=covariance, + ) + return params, emissions, inputs, expected + + +def assert_sample_moments(samples, mean, covariance): + """Check that samples have approximately the expected mean and covariance.""" + num_samples = len(samples) + residuals = np.asarray(samples).reshape(num_samples, -1) - np.asarray(mean).ravel() + whitened = np.linalg.solve(np.linalg.cholesky(covariance), residuals.T).T + np.testing.assert_allclose(whitened.mean(axis=0), 0, atol=5 / np.sqrt(num_samples)) + np.testing.assert_allclose(np.cov(whitened, rowvar=False), np.eye(whitened.shape[1]), + atol=6 / np.sqrt(num_samples)) diff --git a/dynamax/linear_gaussian_ssm/input_timing_test.py b/dynamax/linear_gaussian_ssm/input_timing_test.py new file mode 100644 index 000000000..0a97491cc --- /dev/null +++ b/dynamax/linear_gaussian_ssm/input_timing_test.py @@ -0,0 +1,41 @@ +"""Tests for LGSSM input and parameter indexing and covariance lengths.""" + +import numpy as np +import pytest +from jax import jit, random as jr, vmap + +from dynamax.linear_gaussian_ssm.inference import lgssm_filter, lgssm_posterior_sample, lgssm_smoother +from dynamax.linear_gaussian_ssm.inference_test_utils import make_timing_case, assert_sample_moments + + +@pytest.mark.parametrize("num_timesteps", [1, 3]) +@pytest.mark.parametrize("smoother", [lgssm_smoother]) +def test_inference_input_indexing(num_timesteps, smoother): + """Check input and parameter indexing in LGSSM filtering and smoothing against exact results.""" + params, emissions, inputs, expected = make_timing_case(num_timesteps) + actual = jit(smoother)(params, emissions, inputs)._asdict() + for field, value in actual.items(): + if value is not None: + np.testing.assert_allclose(value, expected[field], atol=2e-5, rtol=2e-5) + + +@pytest.mark.parametrize("num_timesteps", [1, 3]) +@pytest.mark.parametrize("sampler", [lgssm_posterior_sample]) +def test_posterior_sample_input_indexing(num_timesteps, sampler): + """Check input and parameter indexing in LGSSM posterior sampling + against the exact mean and covariance. + """ + params, emissions, inputs, expected = make_timing_case(num_timesteps) + samples = jit(vmap(lambda key: sampler(key, params, emissions, inputs)))( + jr.split(jr.PRNGKey(13), 12000) + ) + assert_sample_moments(samples, expected["smoothed_means"], expected["joint_covariance"]) + + +@pytest.mark.parametrize("filter_fn", [lgssm_filter]) +def test_rejects_short_dynamics_covariance(filter_fn): + """Check that LGSSM filters reject incorrect dynamics covariance lengths.""" + params, emissions, inputs, _ = make_timing_case() + params = params._replace(dynamics=params.dynamics._replace(cov=params.dynamics.cov[:-1])) + with pytest.raises(AssertionError): + filter_fn(params, emissions, inputs) diff --git a/dynamax/linear_gaussian_ssm/models_test.py b/dynamax/linear_gaussian_ssm/models_test.py index 0dc37abc4..ddd2161eb 100644 --- a/dynamax/linear_gaussian_ssm/models_test.py +++ b/dynamax/linear_gaussian_ssm/models_test.py @@ -4,8 +4,9 @@ from functools import partial from itertools import count +import numpy as np import pytest -from jax import vmap +from jax import jit, vmap from jax.tree_util import tree_leaves, tree_map import jax.numpy as jnp import jax.random as jr @@ -22,6 +23,35 @@ ] +@pytest.mark.parametrize("num_forecast_timesteps", [1, 3]) +@pytest.mark.parametrize("has_inputs", [False, True]) +def test_forecast_matches_expected_moments(num_forecast_timesteps, has_inputs): + """Check forecast state and observation means and covariances against hand-computed values.""" + model = LinearGaussianSSM(state_dim=1, emission_dim=1, input_dim=int(has_inputs)) + params, _ = model.initialize( + initial_mean=jnp.zeros(1), initial_covariance=jnp.eye(1), + dynamics_weights=jnp.array([[0.5]]), dynamics_covariance=jnp.eye(1), + dynamics_input_weights=jnp.ones((1, int(has_inputs))), + emission_weights=jnp.eye(1), emission_covariance=jnp.eye(1)) + params = params._replace(dynamics=params.dynamics._replace( + bias=None, input_weights=params.dynamics.input_weights if has_inputs else None)) + inputs = jnp.zeros((1, 1)) if has_inputs else None + future_inputs = jnp.array([[10.0], [20.0], [30.0]])[:num_forecast_timesteps] if has_inputs else None + + forecasts = jit(partial(model.forecast, num_forecast_timesteps=num_forecast_timesteps))( + params, jnp.array([[4.88]]), inputs=inputs, forecast_inputs=future_inputs) + + # Filtering y=4.88 with unit prior/emission variance gives m=2.44, P=0.5. + # Each future step applies m <- 0.5*m + u and P <- 0.25*P + 1. + means = [11.22, 25.61, 42.805] if has_inputs else [1.22, 0.61, 0.305] + expected_means = np.array(means)[:num_forecast_timesteps, None] + expected_covs = np.array([1.125, 1.28125, 1.3203125])[:num_forecast_timesteps, None, None] + np.testing.assert_allclose(forecasts[0], expected_means, rtol=1e-6) + np.testing.assert_allclose(forecasts[1], expected_covs, rtol=1e-6) + np.testing.assert_allclose(forecasts[2], expected_means, rtol=1e-6) + np.testing.assert_allclose(forecasts[3], expected_covs + 1, rtol=1e-6) + + @pytest.mark.parametrize(["cls", "kwargs", "inputs"], CONFIGS) def test_sample_and_fit(cls, kwargs, inputs): """ From b653f40321be24bc36566445afde8419f7916027 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:28:08 -0700 Subject: [PATCH 03/24] fix: correct serial LGSSM input and parameter indexing and forecast start time --- dynamax/linear_gaussian_ssm/inference.py | 39 ++++++++++++++++-------- dynamax/linear_gaussian_ssm/models.py | 16 ++++++---- 2 files changed, 36 insertions(+), 19 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/inference.py b/dynamax/linear_gaussian_ssm/inference.py index 03174d0f1..6e80f2ae0 100644 --- a/dynamax/linear_gaussian_ssm/inference.py +++ b/dynamax/linear_gaussian_ssm/inference.py @@ -36,10 +36,13 @@ class ParamsLGSSMInitial(NamedTuple): class ParamsLGSSMDynamics(NamedTuple): - r"""Parameters of the emission distribution + r"""Parameters of the dynamics distribution - $$p(z_{t+1} \mid z_t, u_t) = \mathcal{N}(z_{t+1} \mid F z_t + B u_t + b, Q)$$ + $$p(z_{t+1} \mid z_t, u_{t+1}) = \mathcal{N}(z_{t+1} \mid F z_t + B u_{t+1} + b, Q)$$ + Time-varying parameter arrays have one entry per timestep. The transition + from z_t to z_{t+1} uses entry t+1. Entry zero is ignored; initial-state + parameters are supplied separately in params.initial. The tuple doubles as a container for the ParameterProperties. :param weights: dynamics weights $F$ @@ -265,7 +268,7 @@ def _condition_on(prior_mean: Float[Array, "state_dim"], emission: Float[Array, "emission_dim"]): r"""Condition a Gaussian potential on a new linear Gaussian observation p(z_t \mid y_t, u_t, y_{1:t-1}, u_{1:t-1}) - propto p(z_t \mid y_{1:t-1}, u_{1:t-1}) p(y_t \mid z_t, u_t) + propto p(z_t \mid y_{1:t-1}, u_{1:t}) p(y_t \mid z_t, u_t) = N(z_t \mid m, P) N(y_t \mid H_t z_t + D_t u_t + d_t, R_t) = N(z_t \mid mm, PP) where @@ -329,6 +332,11 @@ def preprocess_params_and_inputs(params: ParamsLGSSM, assert params.emissions.weights is not None assert params.emissions.cov is not None + for param, ndim in ((params.dynamics.weights, 2), (params.dynamics.input_weights, 2), + (params.dynamics.bias, 1), (params.dynamics.cov, 2)): + if param is not None and not callable(param) and param.ndim == ndim + 1: + assert param.shape[0] == num_timesteps + # Get shapes emission_dim, state_dim = params.emissions.weights.shape[-2:] @@ -486,12 +494,12 @@ def _log_likelihood(pred_mean, pred_cov, H, D, d, R, u, y): return MVNLowRank(m, R, L).log_prob(y) - def _step(carry, t): + def _step(carry, t, predict=True): """Run one step of the Kalman filter.""" ll, pred_mean, pred_cov = carry # Shorthand: get parameters and inputs for time index t - F, B, b, Q, H, D, d, R = _get_params(params, num_timesteps, t) + H, D, d, R = _get_params(params, num_timesteps, t)[4:] u = inputs[t] y = emissions[t] @@ -502,13 +510,18 @@ def _step(carry, t): filtered_mean, filtered_cov = _condition_on(pred_mean, pred_cov, H, D, d, R, u, y) # Predict the next state - pred_mean, pred_cov = _predict(filtered_mean, filtered_cov, F, B, b, Q, u) + if predict: + F, B, b, Q = _get_params(params, num_timesteps, t + 1)[:4] + pred_mean, pred_cov = _predict(filtered_mean, filtered_cov, F, B, b, Q, inputs[t + 1]) return (ll, pred_mean, pred_cov), (filtered_mean, filtered_cov) # Run the Kalman filter carry = (0.0, params.initial.mean, params.initial.cov) - (ll, _, _), (filtered_means, filtered_covs) = lax.scan(_step, carry, jnp.arange(num_timesteps)) + carry, (filtered_means, filtered_covs) = lax.scan(_step, carry, jnp.arange(num_timesteps - 1)) + (ll, _, _), (filtered_mean, filtered_cov) = _step(carry, num_timesteps - 1, predict=False) + filtered_means = jnp.concatenate((filtered_means, filtered_mean[None])) + filtered_covs = jnp.concatenate((filtered_covs, filtered_cov[None])) return PosteriorGSSMFiltered(marginal_loglik=ll, filtered_means=filtered_means, filtered_covariances=filtered_covs) @@ -544,9 +557,9 @@ def _step(carry, args): smoothed_mean_next, smoothed_cov_next = carry t, filtered_mean, filtered_cov = args - # Get parameters and inputs for time index t - F, B, b, Q = _get_params(params, num_timesteps, t)[:4] - u = inputs[t] + # Get parameters and inputs for time index t+1 + F, B, b, Q = _get_params(params, num_timesteps, t + 1)[:4] + u = inputs[t + 1] # This is like the Kalman gain but in reverse # See Eq 8.11 of Saarka's "Bayesian Filtering and Smoothing" @@ -614,9 +627,9 @@ def _step(carry, args): next_state = carry key, filtered_mean, filtered_cov, t = args - # Shorthand: get parameters and inputs for time index t - F, B, b, Q = _get_params(params, num_timesteps, t)[:4] - u = inputs[t] + # Shorthand: get parameters and inputs for time index t+1 + F, B, b, Q = _get_params(params, num_timesteps, t + 1)[:4] + u = inputs[t + 1] # Condition on next state smoothed_mean, smoothed_cov = _condition_on(filtered_mean, filtered_cov, F, B, b, Q, u, next_state) diff --git a/dynamax/linear_gaussian_ssm/models.py b/dynamax/linear_gaussian_ssm/models.py index 332ce3a37..72ce24e85 100644 --- a/dynamax/linear_gaussian_ssm/models.py +++ b/dynamax/linear_gaussian_ssm/models.py @@ -18,6 +18,7 @@ from dynamax.linear_gaussian_ssm.inference import lgssm_joint_sample, lgssm_filter, lgssm_smoother, lgssm_posterior_sample from dynamax.linear_gaussian_ssm.inference import ParamsLGSSM, ParamsLGSSMInitial, ParamsLGSSMDynamics, ParamsLGSSMEmissions from dynamax.linear_gaussian_ssm.inference import PosteriorGSSMFiltered, PosteriorGSSMSmoothed +from dynamax.linear_gaussian_ssm.inference import _predict, preprocess_params_and_inputs from dynamax.parameters import ParameterProperties, ensure_all_or_none_trainable from dynamax.types import PRNGKeyT, Scalar from dynamax.utils.bijectors import RealToPSDBijector @@ -328,13 +329,14 @@ def forecast(self, Float[Array, "num_forecast_timesteps state_dim state_dim"], Float[Array, "num_forecast_timesteps emission_dim"], Float[Array, "num_forecast_timesteps emission_dim emission_dim"]]: - """Compute the marginal filtering distribution for each time step. + """Compute the marginal forecast distribution for each time step. Args: params: model parameters. emissions: sequence of observations. num_forecast_timesteps: number of timesteps to forecast inputs: optional sequence of inputs. + forecast_inputs: optional future inputs, beginning with the first forecast timestep. Returns: forecast state means @@ -346,10 +348,14 @@ def forecast(self, filtered_post = lgssm_filter(params, emissions, inputs) # Forecast into the future from the last timestep + params, forecast_inputs = preprocess_params_and_inputs( + params, num_forecast_timesteps, forecast_inputs) + forecast_mean, forecast_cov = _predict( + filtered_post.filtered_means[-1], filtered_post.filtered_covariances[-1], + params.dynamics.weights, params.dynamics.input_weights, params.dynamics.bias, + params.dynamics.cov, forecast_inputs[0]) forecast_params = ParamsLGSSM( - initial=ParamsLGSSMInitial( - mean=filtered_post.filtered_means[-1], - cov=filtered_post.filtered_covariances[-1]), + initial=ParamsLGSSMInitial(mean=forecast_mean, cov=forecast_cov), dynamics=params.dynamics, emissions=ParamsLGSSMEmissions( weights=params.emissions.weights, @@ -359,8 +365,6 @@ def forecast(self, ) dummy_emissions = jnp.zeros((num_forecast_timesteps, self.emission_dim)) - forecast_inputs = forecast_inputs if forecast_inputs is not None else \ - jnp.zeros((num_forecast_timesteps, 0)) forecast_states = lgssm_filter(forecast_params, dummy_emissions, forecast_inputs) # Forecast future emissions From 987bbaf586018f9ef705c3032b9ccd361d072f7b Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:25:50 -0700 Subject: [PATCH 04/24] test: check LGSSM input indexing in EM and Gibbs sufficient statistics --- dynamax/linear_gaussian_ssm/models_test.py | 68 ++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/dynamax/linear_gaussian_ssm/models_test.py b/dynamax/linear_gaussian_ssm/models_test.py index ddd2161eb..2f9387659 100644 --- a/dynamax/linear_gaussian_ssm/models_test.py +++ b/dynamax/linear_gaussian_ssm/models_test.py @@ -3,15 +3,18 @@ """ from functools import partial from itertools import count +from types import SimpleNamespace import numpy as np import pytest +import jax from jax import jit, vmap from jax.tree_util import tree_leaves, tree_map import jax.numpy as jnp import jax.random as jr from dynamax.linear_gaussian_ssm import LinearGaussianSSM +import dynamax.linear_gaussian_ssm.models as models from dynamax.linear_gaussian_ssm import LinearGaussianConjugateSSM from dynamax.utils.utils import monotonically_increasing @@ -131,3 +134,68 @@ def test_fit_blocked_gibbs_batched(): _, y_obs = vmap(partial(model.sample, params, num_timesteps=num_timesteps))(m_keys) model.fit_blocked_gibbs(next(keys), params, sample_size=6, emissions=y_obs) + + +def test_e_step_input_statistics(monkeypatch): + """Check input indexing in EM dynamics sufficient statistics: pair z_{t-1} with u_t.""" + model = LinearGaussianSSM( + state_dim=1, emission_dim=1, input_dim=1) + inputs = jnp.array([[100.0], [7.0], [-2.0]]) + emissions = jnp.array([[3.0], [-4.0], [8.0]]) + posterior = SimpleNamespace( + smoothed_means=jnp.array([[2.0], [5.0], [-1.0]]), + smoothed_covariances=jnp.array([[[0.5]], [[0.25]], [[0.75]]]), + smoothed_cross_covariances=jnp.array([[[10.1]], [[-5.2]]]), + marginal_loglik=-12.0) + monkeypatch.setattr(models, "lgssm_smoother", lambda *args: posterior) + + (_, dynamics, emission), _ = model.e_step(None, emissions, inputs) + + # Dynamics regressors are [2, 7, 1] and [5, -2, 1], with responses 5 and -1. + # Add posterior variances to squared means. For adjacent-state products, + # use the supplied posterior expectations directly. + expected_xx = np.array([[29.75, 4.0, 7.0], [4.0, 53.0, 5.0], [7.0, 5.0, 2.0]]) + expected_xy = np.array([[4.9], [37.0], [4.0]]) + np.testing.assert_allclose(dynamics[0], expected_xx) + np.testing.assert_allclose(dynamics[1], expected_xy, rtol=1e-6) + np.testing.assert_allclose(dynamics[2], [[27.0]]) + assert dynamics[3] == 2 + + # Emission input statistics still include input[0]. + np.testing.assert_allclose(emission[1][1], [256.0]) + + +def test_gibbs_input_statistics(monkeypatch): + """Check input indexing in Gibbs dynamics sufficient statistics: pair z_{t-1} with u_t.""" + model = LinearGaussianConjugateSSM( + state_dim=1, emission_dim=1, input_dim=1) + params, _ = model.initialize(jr.PRNGKey(0)) + emissions = jnp.array([[[2.0], [5.0], [-1.0]], [[-3.0], [4.0], [6.0]]]) + inputs = jnp.array([[[100.0], [7.0], [-2.0]], [[200.0], [-5.0], [3.0]]]) + # Fix the latent-state samples to check input indexing in Gibbs sufficient + # statistics independently of posterior sampling. + monkeypatch.setattr(models, "lgssm_posterior_sample", lambda key, params, emissions, inputs: emissions) + updates = [] + + def record_update(prior, stats): + """Record the sufficient statistics passed to the parameter update.""" + updates.append(stats) + # Only the statistics are tested; return the prior so sampling can continue. + return prior + + monkeypatch.setattr(models, "mniw_posterior_update", record_update) + with jax.disable_jit(): + model.fit_blocked_gibbs(jr.PRNGKey(1), params, sample_size=1, emissions=emissions, inputs=inputs) + + # Four transitions, with regressors [2,7,1], [5,-2,1], [-3,-5,1], [4,3,1]. + # The inputs 100 and 200 affect no transition and cannot enter these products. + expected_xx = np.array([[54.0, 31.0, 8.0], [31.0, 87.0, 3.0], [8.0, 3.0, 4.0]]) + expected_xy = np.array([[17.0], [35.0], [14.0]]) + dynamics = updates[0] + np.testing.assert_allclose(dynamics[0], expected_xx) + np.testing.assert_allclose(dynamics[1], expected_xy) + np.testing.assert_allclose(dynamics[2], [[78.0]]) + assert dynamics[3] == 4 + + # Emission input statistics still include each sequence's first input. + np.testing.assert_allclose(updates[1][1][1], [-365.0]) From f2bb78c8974eb1faba21e752c38c1bbeabb77b64 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:27:22 -0700 Subject: [PATCH 05/24] fix: correct LGSSM input indexing in EM sufficient statistics --- dynamax/linear_gaussian_ssm/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/models.py b/dynamax/linear_gaussian_ssm/models.py index 72ce24e85..643793648 100644 --- a/dynamax/linear_gaussian_ssm/models.py +++ b/dynamax/linear_gaussian_ssm/models.py @@ -415,7 +415,7 @@ def e_step(self, # Append bias to the inputs inputs = jnp.concatenate((inputs, jnp.ones((num_timesteps, 1))), axis=1) - up = inputs[:-1] + up = inputs[1:] u = inputs y = emissions @@ -425,7 +425,7 @@ def e_step(self, init_stats = (Ex0, Ex0x0T, 1) # expected sufficient statistics for the dynamics tfd.Distribution - # let zp[t] = [x[t], u[t]] for t = 0...T-2 + # let zp[t] = [x[t], u[t+1]] for t = 0...T-2 # let xn[t] = x[t+1] for t = 0...T-2 sum_zpzpT = jnp.block([[Exp.T @ Exp, Exp.T @ up], [up.T @ Exp, up.T @ up]]) sum_zpzpT = sum_zpzpT.at[:self.state_dim, :self.state_dim].add(Vxp.sum(0)) From d3d0de4a6f94cd73aca7b84cff50a8ffa53b3c84 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:27:23 -0700 Subject: [PATCH 06/24] fix: correct LGSSM input indexing in Gibbs sufficient statistics --- dynamax/linear_gaussian_ssm/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/models.py b/dynamax/linear_gaussian_ssm/models.py index 643793648..cfebeffd8 100644 --- a/dynamax/linear_gaussian_ssm/models.py +++ b/dynamax/linear_gaussian_ssm/models.py @@ -674,12 +674,12 @@ def sufficient_stats_from_sample(y, inputs, states): inputs_joint = jnp.concatenate((inputs, jnp.ones((num_timesteps, 1))), axis=1) # Let xn[t] = x[t+1] for t = 0...T-2 x, xp, xn = states, states[:-1], states[1:] - u, up = inputs_joint, inputs_joint[:-1] + u, up = inputs_joint, inputs_joint[1:] init_stats = (x[0], jnp.outer(x[0], x[0]), 1) # Quantities for the dynamics distribution - # Let zp[t] = [x[t], u[t]] for t = 0...T-2 + # Let zp[t] = [x[t], u[t+1]] for t = 0...T-2 sum_zpzpT = jnp.block([[xp.T @ xp, xp.T @ up], [up.T @ xp, up.T @ up]]) sum_zpxnT = jnp.block([[xp.T @ xn], [up.T @ xn]]) sum_xnxnT = xn.T @ xn From 0d3d39b9c2aad993d0e81308b1dde5cd0a475d8c Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:25:50 -0700 Subject: [PATCH 07/24] test: check parallel and information-form LGSSM input and parameter indexing and covariance lengths --- .../info_inference_test.py | 22 ++++++++++++- .../linear_gaussian_ssm/input_timing_test.py | 33 +++++++++++++++++-- 2 files changed, 51 insertions(+), 4 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/info_inference_test.py b/dynamax/linear_gaussian_ssm/info_inference_test.py index d9ab8672c..d3ed6677d 100644 --- a/dynamax/linear_gaussian_ssm/info_inference_test.py +++ b/dynamax/linear_gaussian_ssm/info_inference_test.py @@ -4,12 +4,13 @@ import jax.numpy as jnp from functools import partial +from types import SimpleNamespace from jax import random as jr from dynamax.linear_gaussian_ssm.inference import lgssm_smoother, lgssm_filter from dynamax.linear_gaussian_ssm.inference import ParamsLGSSM, ParamsLGSSMInitial, ParamsLGSSMDynamics, ParamsLGSSMEmissions from dynamax.linear_gaussian_ssm.info_inference import lgssm_info_filter, lgssm_info_smoother, info_to_moment_form -from dynamax.linear_gaussian_ssm.info_inference import ParamsLGSSMInfo +from dynamax.linear_gaussian_ssm.info_inference import ParamsLGSSMInfo, lds_to_block_tridiag from dynamax.utils.utils import has_tpu # Use lower tolerance for TPU tests. @@ -180,3 +181,22 @@ def test_filtered_means(self): def test_filtered_covs(self): """Test filtered covariances.""" assert allclose(self.info_filtered_covs, self.lgssm_moment_posterior.filtered_covariances) + + +def test_block_tridiag_transition_indices(): + """Check that transitions from z_{t-1} to z_t use parameters and inputs at t.""" + weights = jnp.array([jnp.nan, 2., 3.])[:, None, None] + input_weights = jnp.array([jnp.nan, 5., 7.])[:, None, None] + covariances = jnp.array([jnp.nan, 11., 13.])[:, None, None] + lds = SimpleNamespace( + initial_mean=jnp.zeros(1), initial_covariance=jnp.eye(1), + dynamics_matrix=weights.__getitem__, dynamics_input_weights=input_weights.__getitem__, + dynamics_noise_covariance=covariances.__getitem__, + emissions_matrix=lambda t: jnp.eye(1), emissions_input_weights=lambda t: jnp.zeros((1, 1)), + emissions_noise_covariance=lambda t: jnp.eye(1), + ) + diagonal, lower, h = lds_to_block_tridiag(lds, jnp.zeros((3, 1)), jnp.array([[100.], [2.], [-1.]])) + # Expand z0**2 + (z1-2*z0-10)**2/11 + (z2-3*z1+7)**2/13 + sum(zt**2). + assert allclose(diagonal[:, 0, 0], jnp.array([2 + 4/11, 1 + 1/11 + 9/13, 1 + 1/13])) + assert allclose(lower[:, 0, 0], jnp.array([-2/11, -3/13])) + assert allclose(h[:, 0], jnp.array([-20/11, 10/11 + 21/13, -7/13])) diff --git a/dynamax/linear_gaussian_ssm/input_timing_test.py b/dynamax/linear_gaussian_ssm/input_timing_test.py index 0a97491cc..41b634791 100644 --- a/dynamax/linear_gaussian_ssm/input_timing_test.py +++ b/dynamax/linear_gaussian_ssm/input_timing_test.py @@ -2,25 +2,50 @@ import numpy as np import pytest +import jax.numpy as jnp from jax import jit, random as jr, vmap from dynamax.linear_gaussian_ssm.inference import lgssm_filter, lgssm_posterior_sample, lgssm_smoother +from dynamax.linear_gaussian_ssm import parallel_lgssm_filter, parallel_lgssm_posterior_sample, parallel_lgssm_smoother +from dynamax.linear_gaussian_ssm.info_inference import ParamsLGSSMInfo, info_to_moment_form, lgssm_info_filter, lgssm_info_smoother from dynamax.linear_gaussian_ssm.inference_test_utils import make_timing_case, assert_sample_moments +def _to_info_form(params): + """Express the same time-varying model using precisions instead of covariances.""" + return ParamsLGSSMInfo( + initial_mean=params.initial.mean, + initial_precision=jnp.linalg.inv(params.initial.cov), + dynamics_weights=params.dynamics.weights, + dynamics_precision=jnp.linalg.inv(params.dynamics.cov), + dynamics_input_weights=params.dynamics.input_weights, + dynamics_bias=params.dynamics.bias, + emission_weights=params.emissions.weights, + emission_precision=jnp.linalg.inv(params.emissions.cov), + emission_input_weights=params.emissions.input_weights, + emission_bias=params.emissions.bias, + ) + + @pytest.mark.parametrize("num_timesteps", [1, 3]) -@pytest.mark.parametrize("smoother", [lgssm_smoother]) +@pytest.mark.parametrize("smoother", [lgssm_smoother, parallel_lgssm_smoother, lgssm_info_smoother]) def test_inference_input_indexing(num_timesteps, smoother): """Check input and parameter indexing in LGSSM filtering and smoothing against exact results.""" params, emissions, inputs, expected = make_timing_case(num_timesteps) + if smoother == lgssm_info_smoother: + params = _to_info_form(params) actual = jit(smoother)(params, emissions, inputs)._asdict() + if smoother == lgssm_info_smoother: + for field in ("filtered", "smoothed"): + actual[f"{field}_means"], actual[f"{field}_covariances"] = info_to_moment_form( + actual.pop(f"{field}_etas"), actual.pop(f"{field}_precisions")) for field, value in actual.items(): if value is not None: np.testing.assert_allclose(value, expected[field], atol=2e-5, rtol=2e-5) @pytest.mark.parametrize("num_timesteps", [1, 3]) -@pytest.mark.parametrize("sampler", [lgssm_posterior_sample]) +@pytest.mark.parametrize("sampler", [lgssm_posterior_sample, parallel_lgssm_posterior_sample]) def test_posterior_sample_input_indexing(num_timesteps, sampler): """Check input and parameter indexing in LGSSM posterior sampling against the exact mean and covariance. @@ -32,10 +57,12 @@ def test_posterior_sample_input_indexing(num_timesteps, sampler): assert_sample_moments(samples, expected["smoothed_means"], expected["joint_covariance"]) -@pytest.mark.parametrize("filter_fn", [lgssm_filter]) +@pytest.mark.parametrize("filter_fn", [lgssm_filter, parallel_lgssm_filter, lgssm_info_filter]) def test_rejects_short_dynamics_covariance(filter_fn): """Check that LGSSM filters reject incorrect dynamics covariance lengths.""" params, emissions, inputs, _ = make_timing_case() params = params._replace(dynamics=params.dynamics._replace(cov=params.dynamics.cov[:-1])) + if filter_fn == lgssm_info_filter: + params = _to_info_form(params) with pytest.raises(AssertionError): filter_fn(params, emissions, inputs) From 9fb889445c29601053e1957787d04eef5feaefda Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:28:08 -0700 Subject: [PATCH 08/24] fix: correct input and parameter indexing in parallel LGSSM filtering and smoothing --- .../linear_gaussian_ssm/parallel_inference.py | 26 +++++++++++-------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/parallel_inference.py b/dynamax/linear_gaussian_ssm/parallel_inference.py index e9cba6d19..bda774288 100644 --- a/dynamax/linear_gaussian_ssm/parallel_inference.py +++ b/dynamax/linear_gaussian_ssm/parallel_inference.py @@ -159,6 +159,10 @@ def _initialize_filtering_messages( """Preprocess observations to construct input for filtering assocative scan.""" num_timesteps = emissions.shape[0] + for param, ndim in ((params.dynamics.weights, 2), (params.dynamics.input_weights, 2), + (params.dynamics.bias, 1), (params.dynamics.cov, 2)): + if param is not None and not callable(param) and param.ndim == ndim + 1: + assert param.shape[0] == num_timesteps inputs = _zeros_if_none(inputs, (num_timesteps, 0)) # Get the emission covariance dimension @@ -187,13 +191,13 @@ def _first_message(params, y, u): return A, b, C, J, eta, logZ - @partial(vmap, in_axes=(None, 0, 0, 0, 0)) - def _generic_message(params, y, ut, utm1, t): + @partial(vmap, in_axes=(None, 0, 0, 0)) + def _generic_message(params, y, ut, t): """Compute the generic filtering message.""" - F = _get_one_param(params.dynamics.weights, 2, t-1) - B = _get_one_param(params.dynamics.input_weights, 2, t-1) - b = _get_one_param(params.dynamics.bias, 1, t-1) - Q = _get_one_param(params.dynamics.cov, 2, t-1) + F = _get_one_param(params.dynamics.weights, 2, t) + B = _get_one_param(params.dynamics.input_weights, 2, t) + b = _get_one_param(params.dynamics.bias, 1, t) + Q = _get_one_param(params.dynamics.cov, 2, t) H = _get_one_param(params.emissions.weights, 2, t) D = _get_one_param(params.emissions.input_weights, 2, t) d = _get_one_param(params.emissions.bias, 1, t) @@ -202,10 +206,10 @@ def _generic_message(params, y, ut, utm1, t): S_inv = _emissions_scale(Q, H, R) K = Q @ H.T @ S_inv - bias_tm1 = B @ utm1 + b - innov = (y - D @ ut - d - H @ bias_tm1) + bias_t = B @ ut + b + innov = (y - D @ ut - d - H @ bias_t) A = F - K @ H @ F - b = bias_tm1 + K @ innov + b = bias_t + K @ innov C = symmetrize(Q - K @ H @ Q) eta = F.T @ H.T @ S_inv @ innov J = symmetrize(F.T @ H.T @ S_inv @ H @ F) @@ -214,7 +218,7 @@ def _generic_message(params, y, ut, utm1, t): return A, b, C, J, eta, logZ A0, b0, C0, J0, eta0, logZ0 = _first_message(params, emissions[0], inputs[0]) - At, bt, Ct, Jt, etat, logZt = _generic_message(params, emissions[1:], inputs[1:], inputs[:-1], jnp.arange(1, len(emissions))) + At, bt, Ct, Jt, etat, logZt = _generic_message(params, emissions[1:], inputs[1:], jnp.arange(1, len(emissions))) return FilterMessage( A=jnp.concatenate([A0[None], At]), @@ -314,7 +318,7 @@ def _generic_message(params, m, P, ut, t): return E, g, L En, gn, Ln = _last_message(filtered_means[-1], filtered_covariances[-1]) - Et, gt, Lt = _generic_message(params, filtered_means[:-1], filtered_covariances[:-1], inputs[:-1], jnp.arange(len(filtered_means)-1)) + Et, gt, Lt = _generic_message(params, filtered_means[:-1], filtered_covariances[:-1], inputs[1:], jnp.arange(1, len(filtered_means))) return SmoothMessage( E=jnp.concatenate([Et, En[None]]), From de81fdec18982bbfc9d398e90c1765fd285dc29a Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:28:08 -0700 Subject: [PATCH 09/24] fix: pass missing inputs to parallel LGSSM backward sampling --- dynamax/linear_gaussian_ssm/parallel_inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dynamax/linear_gaussian_ssm/parallel_inference.py b/dynamax/linear_gaussian_ssm/parallel_inference.py index bda774288..d7bd65a71 100644 --- a/dynamax/linear_gaussian_ssm/parallel_inference.py +++ b/dynamax/linear_gaussian_ssm/parallel_inference.py @@ -385,7 +385,7 @@ def _initialize_sampling_messages(key, params, filtered_means, filtered_covarian Given parallel smoothing messages `z_i ~ N(E_i z_{i+1} + g_i, L_i)`, the parallel sampling messages are `(E_i,h_i)` where `h_i ~ N(g_i, L_i)`. """ - E, g, L = _initialize_smoothing_messages(params, filtered_means, filtered_covariances) + E, g, L = _initialize_smoothing_messages(params, filtered_means, filtered_covariances, inputs) return SampleMessage(E=E, h=MVN(g, L).sample(seed=key)) From 7bc156c9e9914e343cec032d1d0592b50fbc3b86 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:28:08 -0700 Subject: [PATCH 10/24] fix: correct input and parameter indexing in information-form LGSSM inference --- dynamax/linear_gaussian_ssm/info_inference.py | 54 +++++++++++-------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/info_inference.py b/dynamax/linear_gaussian_ssm/info_inference.py index fe62afd8e..7d31d7328 100644 --- a/dynamax/linear_gaussian_ssm/info_inference.py +++ b/dynamax/linear_gaussian_ssm/info_inference.py @@ -96,16 +96,16 @@ def _info_predict(eta, Lambda, F, Q_prec, B, u, b): Marginalising over the uncertainty in z_t the predicted latent state at the next time step is given by: - p(z_{t+1}\mid z_t, u_t) - = \int p(z_{t+1}, z_t \mid u_t) dz_t - = \int N(z_t \mid mu_t, Sigma_t) N(z_{t+1} \mid F z_t + B u_t + b, Q) dz_t + p(z_{t+1}\mid z_t, u_{t+1}) + = \int p(z_{t+1}, z_t \mid u_{t+1}) dz_t + = \int N(z_t \mid mu_t, Sigma_t) N(z_{t+1} \mid F z_t + B u_{t+1} + b, Q) dz_t = N(z_t \mid m_{t+1\midt}, Sigma_{t+1\mid t}) with - m_{t+1 \mid t} = F m_t + B u_t + b + m_{t+1 \mid t} = F m_t + B u_{t+1} + b Sigma_{t+1 \mid t} = F Sigma_t F^T + Q The corresponding information form parameters are: - eta_{t+1 \mid t} = K eta_t + Lambda_{t+1 \mid t} (B u_t + b) + eta_{t+1 \mid t} = K eta_t + Lambda_{t+1 \mid t} (B u_{t+1} + b) Lambda_{t+1 \mid t} = L Q_prec L^T + K Lambda_t K^T where K = Q_prec F ( Lambda_t + F^T Q_prec F)^{-1} @@ -187,18 +187,18 @@ def lgssm_info_filter( filtered_precisions """ num_timesteps = len(emissions) + for param, ndim in ((params.dynamics_weights, 2), (params.dynamics_input_weights, 2), + (params.dynamics_bias, 1), (params.dynamics_precision, 2)): + if param is not None and not callable(param) and param.ndim == ndim + 1: + assert param.shape[0] == num_timesteps inputs = jnp.zeros((num_timesteps, 0)) if inputs is None else inputs - def _filter_step(carry, t): + def _filter_step(carry, t, predict=True): """Run a single step of the Kalman filter.""" ll, pred_eta, pred_prec = carry # Shorthand: get parameters and inputs for time index t - F = _get_params(params.dynamics_weights, 2, t) - Q_prec = _get_params(params.dynamics_precision, 2, t) H = _get_params(params.emission_weights, 2, t) R_prec = _get_params(params.emission_precision, 2, t) - B = _get_params(params.dynamics_input_weights, 2, t) - b = _get_params(params.dynamics_bias, 1, t) D = _get_params(params.emission_input_weights, 2, t) d = _get_params(params.emission_bias, 1, t) u = inputs[t] @@ -212,14 +212,22 @@ def _filter_step(carry, t): filtered_eta, filtered_prec = _info_condition_on(pred_eta, pred_prec, H, R_prec, D, u, d, y) # Predict the next state - pred_eta, pred_prec = _info_predict(filtered_eta, filtered_prec, F, Q_prec, B, u, b) + if predict: + F = _get_params(params.dynamics_weights, 2, t + 1) + Q_prec = _get_params(params.dynamics_precision, 2, t + 1) + B = _get_params(params.dynamics_input_weights, 2, t + 1) + b = _get_params(params.dynamics_bias, 1, t + 1) + pred_eta, pred_prec = _info_predict(filtered_eta, filtered_prec, F, Q_prec, B, inputs[t + 1], b) return (ll, pred_eta, pred_prec), (filtered_eta, filtered_prec) # Run the Kalman filter initial_eta = params.initial_precision @ params.initial_mean carry = (0.0, initial_eta, params.initial_precision) - (ll, _, _), (filtered_etas, filtered_precisions) = lax.scan(_filter_step, carry, jnp.arange(num_timesteps)) + carry, (filtered_etas, filtered_precisions) = lax.scan(_filter_step, carry, jnp.arange(num_timesteps - 1)) + (ll, _, _), (filtered_eta, filtered_prec) = _filter_step(carry, num_timesteps - 1, predict=False) + filtered_etas = jnp.concatenate((filtered_etas, filtered_eta[None])) + filtered_precisions = jnp.concatenate((filtered_precisions, filtered_prec[None])) return PosteriorGSSMInfoFiltered(marginal_loglik=ll, filtered_etas=filtered_etas, filtered_precisions=filtered_precisions) @@ -255,12 +263,12 @@ def _smooth_step(carry, args): smoothed_eta_next, smoothed_prec_next = carry t, filtered_eta, filtered_prec = args - # Shorthand: get parameters and inputs for time index t - F = _get_params(params.dynamics_weights, 2, t) - B = _get_params(params.dynamics_input_weights, 2, t) - b = _get_params(params.dynamics_bias, 1, t) - Q_prec = _get_params(params.dynamics_precision, 2, t) - u = inputs[t] + # Shorthand: get parameters and inputs for time index t+1 + F = _get_params(params.dynamics_weights, 2, t + 1) + B = _get_params(params.dynamics_input_weights, 2, t + 1) + b = _get_params(params.dynamics_bias, 1, t + 1) + Q_prec = _get_params(params.dynamics_precision, 2, t + 1) + u = inputs[t + 1] # Predict the next state # TODO: Pass predicted params from lgssm_info_filter? @@ -464,16 +472,16 @@ def lds_to_block_tridiag(lds: ParamsLGSSM, # diagonal blocks of precision matrix J_diag = jnp.array([jnp.dot(C(t).T, psd_solve(R(t), C(t))) for t in range(T)]) J_diag = J_diag.at[0].add(jnp.linalg.inv(Q0)) - J_diag = J_diag.at[:-1].add(jnp.array([jnp.dot(A(t).T, psd_solve(Q(t), A(t))) for t in range(T - 1)])) - J_diag = J_diag.at[1:].add(jnp.array([jnp.linalg.inv(Q(t)) for t in range(0, T - 1)])) + J_diag = J_diag.at[:-1].add(jnp.array([jnp.dot(A(t).T, psd_solve(Q(t), A(t))) for t in range(1, T)])) + J_diag = J_diag.at[1:].add(jnp.array([jnp.linalg.inv(Q(t)) for t in range(1, T)])) # lower diagonal blocks of precision matrix - J_lower_diag = jnp.array([-psd_solve(Q(t), A(t)) for t in range(T - 1)]) + J_lower_diag = jnp.array([-psd_solve(Q(t), A(t)) for t in range(1, T)]) # linear potential h = jnp.array([jnp.dot(data[t] - D(t) @ inputs[t], psd_solve(R(t), C(t))) for t in range(T)]) h = h.at[0].add(psd_solve(Q0, m0)) - h = h.at[:-1].add(jnp.array([-jnp.dot(A(t).T, psd_solve(Q(t), B(t) @ inputs[t])) for t in range(T - 1)])) - h = h.at[1:].add(jnp.array([psd_solve(Q(t), B(t) @ inputs[t]) for t in range(T - 1)])) + h = h.at[:-1].add(jnp.array([-jnp.dot(A(t).T, psd_solve(Q(t), B(t) @ inputs[t])) for t in range(1, T)])) + h = h.at[1:].add(jnp.array([psd_solve(Q(t), B(t) @ inputs[t]) for t in range(1, T)])) return J_diag, J_lower_diag, h From 64d48683bf06ee1669be1442c8d72a6d804d165b Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:28:09 -0700 Subject: [PATCH 11/24] test: reuse the exact LGSSM example for EKF, UKF, and GGSSM tests --- .../nonlinear_gaussian_ssm/inference_test_utils.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/dynamax/nonlinear_gaussian_ssm/inference_test_utils.py b/dynamax/nonlinear_gaussian_ssm/inference_test_utils.py index 3cf0b879a..ca9c3c967 100644 --- a/dynamax/nonlinear_gaussian_ssm/inference_test_utils.py +++ b/dynamax/nonlinear_gaussian_ssm/inference_test_utils.py @@ -10,6 +10,7 @@ import tensorflow_probability.substrates.jax as tfp from dynamax.linear_gaussian_ssm import LinearGaussianSSM +from dynamax.linear_gaussian_ssm.inference_test_utils import make_timing_case from dynamax.linear_gaussian_ssm.inference import ParamsLGSSM, ParamsLGSSMInitial, ParamsLGSSMDynamics, ParamsLGSSMEmissions from dynamax.nonlinear_gaussian_ssm.models import ParamsNLGSSM, NonlinearGaussianSSM from dynamax.parameters import ParameterProperties @@ -186,3 +187,16 @@ def random_nlgssm_args(key=0, num_timesteps=15, state_dim=4, emission_dim=2): states, emissions = model.sample(params, sample_key, num_timesteps) return params, states, emissions + +def make_affine_timing_case(num_timesteps=3): + """Convert the model from make_timing_case to ParamsNLGSSM.""" + params, emissions, inputs, expected = make_timing_case(num_timesteps) + nonlinear_params = ParamsNLGSSM( + initial_mean=params.initial.mean, + initial_covariance=params.initial.cov, + dynamics_function=lambda state, inpt: inpt * state / 2 + inpt * (inpt + 1) / 2, + dynamics_covariance=params.dynamics.cov, + emission_function=lambda state, inpt: state + inpt / 2, + emission_covariance=params.emissions.cov, + ) + return nonlinear_params, emissions, inputs, expected From 757df05a94cb18fca84daf7be8f4e3159e31e0a4 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 12/24] test: check EKF input and parameter indexing, prediction lengths, and covariance lengths --- ...nference_ekf_test.py => inference_test.py} | 94 ++++++++++++++++++- 1 file changed, 91 insertions(+), 3 deletions(-) rename dynamax/nonlinear_gaussian_ssm/{inference_ekf_test.py => inference_test.py} (53%) diff --git a/dynamax/nonlinear_gaussian_ssm/inference_ekf_test.py b/dynamax/nonlinear_gaussian_ssm/inference_test.py similarity index 53% rename from dynamax/nonlinear_gaussian_ssm/inference_ekf_test.py rename to dynamax/nonlinear_gaussian_ssm/inference_test.py index 3fbf8fccf..695a5ed23 100644 --- a/dynamax/nonlinear_gaussian_ssm/inference_ekf_test.py +++ b/dynamax/nonlinear_gaussian_ssm/inference_test.py @@ -1,16 +1,20 @@ """ Tests for the extended Kalman filter and smoother. """ +import numpy as np +import pytest import jax.numpy as jnp import jax.random as jr from functools import partial from dynamax.linear_gaussian_ssm import lgssm_filter, lgssm_smoother, lgssm_posterior_sample +from dynamax.linear_gaussian_ssm.inference_test_utils import assert_sample_moments from dynamax.nonlinear_gaussian_ssm.inference_ekf import extended_kalman_filter, extended_kalman_smoother, extended_kalman_posterior_sample -from dynamax.nonlinear_gaussian_ssm.inference_test_utils import lgssm_to_nlgssm, random_lgssm_args, random_nlgssm_args +from dynamax.nonlinear_gaussian_ssm.inference_test_utils import lgssm_to_nlgssm, random_lgssm_args, random_nlgssm_args, make_affine_timing_case +from dynamax.nonlinear_gaussian_ssm.models import ParamsNLGSSM from dynamax.nonlinear_gaussian_ssm.sarkka_lib import ekf, eks from dynamax.utils.utils import has_tpu -from jax import vmap +from jax import jit, vmap if has_tpu(): # TPU has very poor numerical stability @@ -128,4 +132,88 @@ def test_extended_kalman_sampler_nonlinear(key=0, num_timesteps=15, sample_size= threshold = 6 * jnp.sqrt(post_variance / sample_size) empirical_means = ekf_samples.mean(0) assert jnp.all(abs(empirical_means - ekf_post.smoothed_means) < threshold) - \ No newline at end of file + + +@pytest.mark.parametrize( + "filter_fn,smoother_fn", + [ + (extended_kalman_filter, extended_kalman_smoother), + ], + ids=["extended"], +) +@pytest.mark.parametrize("num_timesteps", [1, 3]) +def test_kalman_inference_input_indexing(filter_fn, smoother_fn, num_timesteps): + """Check EKF input and parameter indexing by comparing filtering, smoothing, + and prediction results with exact LGSSM results. + """ + params, emissions, inputs, expected = make_affine_timing_case(num_timesteps) + filtered = jit(lambda y, u: filter_fn(params, y, inputs=u))(emissions, inputs) + smoothed = jit(lambda y, u: smoother_fn(params, y, inputs=u))(emissions, inputs) + for posterior in (filtered, smoothed): + for field in posterior._fields: + value = getattr(posterior, field) + if value is not None: + np.testing.assert_allclose(value, expected[field], atol=2e-5, rtol=2e-5) + + +@pytest.mark.parametrize( + "filter_fn", + [extended_kalman_filter], + ids=["extended"], +) +@pytest.mark.parametrize("num_timesteps", [1, 4]) +@pytest.mark.parametrize("time_varying_covariance", [False, True]) +def test_kalman_prediction_length_without_inputs(filter_fn, num_timesteps, time_varying_covariance): + """Check EKF prediction lengths and values without inputs, for static and + time-varying dynamics covariance. + + Check both full-output and prediction-only requests. + """ + covariance = jnp.eye(1) + if time_varying_covariance: + covariance = jnp.broadcast_to(covariance, (num_timesteps, 1, 1)).at[0].set(jnp.nan) + params = ParamsNLGSSM( + initial_mean=jnp.array([2.0]), initial_covariance=jnp.eye(1), + dynamics_function=lambda state: 0.5 * state, dynamics_covariance=covariance, + emission_function=lambda state: state, emission_covariance=jnp.eye(1), + ) + emissions = jnp.array([[0.0], [1.0], [-1.0], [2.0]])[:num_timesteps] + num_predictions = num_timesteps - int(time_varying_covariance) + # Scalar Kalman updates followed by m <- m/2, P <- P/4 + 1. + expected_means = np.array([1/2, 13/34, -51/290, 606/1237])[:num_predictions, None] + expected_covs = np.array([9/8, 77/68, 657/580, 5605/4948])[:num_predictions, None, None] + for output_fields in ( + ["filtered_means", "filtered_covariances", "predicted_means", "predicted_covariances"], + ["predicted_means", "predicted_covariances"], + ): + posterior = jit(lambda y: filter_fn(params, y, output_fields=output_fields))(emissions) + if "filtered_means" not in output_fields: + assert posterior.filtered_means is None + assert posterior.filtered_covariances is None + np.testing.assert_allclose(posterior.predicted_means, expected_means, atol=2e-5, rtol=2e-5) + np.testing.assert_allclose(posterior.predicted_covariances, expected_covs, atol=2e-5, rtol=2e-5) + + +@pytest.mark.parametrize("num_timesteps", [1, 3]) +def test_extended_kalman_posterior_sample_input_indexing(num_timesteps): + """Check input and parameter indexing in EKF posterior sampling + against the exact LGSSM mean and covariance. + """ + params, emissions, inputs, expected = make_affine_timing_case(num_timesteps) + samples = jit(vmap(lambda key: extended_kalman_posterior_sample(key, params, emissions, inputs)))( + jr.split(jr.PRNGKey(31), 10000) + ) + assert_sample_moments(samples, expected["smoothed_means"], expected["joint_covariance"]) + + +def test_kalman_inference_rejects_short_dynamics_covariance(): + """Check that the EKF filter and smoother reject incorrect dynamics covariance lengths.""" + params, emissions, inputs, _ = make_affine_timing_case() + filtered = extended_kalman_filter(params, emissions, inputs=inputs) + params = params._replace(dynamics_covariance=params.dynamics_covariance[:-1]) + for inference in ( + lambda y, u: extended_kalman_filter(params, y, inputs=u), + lambda y, u: extended_kalman_smoother(params, y, filtered_posterior=filtered, inputs=u), + ): + with pytest.raises(AssertionError): + inference(emissions, inputs) From acfe576c768ee37bee9696006a08b660fd817b51 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 13/24] fix: correct EKF input and parameter indexing and prediction lengths --- .../nonlinear_gaussian_ssm/inference_ekf.py | 35 ++++++++++++------- 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/dynamax/nonlinear_gaussian_ssm/inference_ekf.py b/dynamax/nonlinear_gaussian_ssm/inference_ekf.py index 4e815b163..71a553060 100644 --- a/dynamax/nonlinear_gaussian_ssm/inference_ekf.py +++ b/dynamax/nonlinear_gaussian_ssm/inference_ekf.py @@ -53,7 +53,7 @@ def _condition_on(prior_mean: Float[Array, " state_dim"], r"""Condition a Gaussian potential on a new observation. p(z_t | y_t, u_t, y_{1:t-1}, u_{1:t-1}) - propto p(z_t | y_{1:t-1}, u_{1:t-1}) p(y_t | z_t, u_t) + propto p(z_t | y_{1:t-1}, u_{1:t}) p(y_t | z_t, u_t) = N(z_t | m, S) N(y_t | h_t(z_t, u_t), R_t) = N(z_t | mm, SS) where @@ -106,10 +106,14 @@ def extended_kalman_filter(params: ParamsNLGSSM, "predicted_means", "predicted_covariances", and "marginal_loglik". Returns: - post: posterior object. + post: posterior object. Predicted entry ``t`` is for state ``t + 1``. + Predicted fields have length ``T`` when inputs are omitted and + dynamics covariance is static, otherwise ``T - 1``. """ num_timesteps = len(emissions) + predict_final = inputs is None and params.dynamics_covariance.ndim == 2 + assert params.dynamics_covariance.ndim != 3 or params.dynamics_covariance.shape[0] == num_timesteps # Dynamics and emission functions and their Jacobians f, h = params.dynamics_function, params.emission_function @@ -117,12 +121,11 @@ def extended_kalman_filter(params: ParamsNLGSSM, f, h, F, H = (_process_fn(fn, inputs) for fn in (f, h, F, H)) inputs = _process_input(inputs, num_timesteps) - def _step(carry, t): + def _step(carry, t, predict=True): """Iteratively update the state estimate and log likelihood.""" ll, pred_mean, pred_cov = carry # Get parameters and inputs for time index t - Q = _get_params(params.dynamics_covariance, 2, t) R = _get_params(params.emission_covariance, 2, t) u = inputs[t] y = emissions[t] @@ -135,7 +138,10 @@ def _step(carry, t): filtered_mean, filtered_cov = _condition_on(pred_mean, pred_cov, h, H, R, u, y, num_iter) # Predict the next state - pred_mean, pred_cov = _predict(filtered_mean, filtered_cov, f, F, Q, u) + if predict: + Q = _get_params(params.dynamics_covariance, 2, t + 1) + u_next = inputs[t + 1] if not predict_final else u + pred_mean, pred_cov = _predict(filtered_mean, filtered_cov, f, F, Q, u_next) # Build carry and output states carry = (ll, pred_mean, pred_cov) @@ -152,7 +158,11 @@ def _step(carry, t): # Run the extended Kalman filter carry = (0.0, params.initial_mean, params.initial_covariance) - (ll, _, _), outputs = lax.scan(_step, carry, jnp.arange(num_timesteps)) + carry, outputs = lax.scan(_step, carry, jnp.arange(num_timesteps - 1)) + (ll, _, _), final_outputs = _step(carry, num_timesteps - 1, predict=predict_final) + outputs = {key: jnp.concatenate((val, final_outputs[key][None])) for key, val in outputs.items()} + if not predict_final: + outputs = {key: val[:-1] if key.startswith("predicted") else val for key, val in outputs.items()} outputs = {"marginal_loglik": ll, **outputs} posterior_filtered = PosteriorGSSMFiltered( **outputs, @@ -178,6 +188,7 @@ def extended_kalman_smoother(params: ParamsNLGSSM, """ num_timesteps = len(emissions) + assert params.dynamics_covariance.ndim != 3 or params.dynamics_covariance.shape[0] == num_timesteps # Get filtered posterior if filtered_posterior is None: @@ -198,10 +209,10 @@ def _step(carry, args): smoothed_mean_next, smoothed_cov_next = carry t, filtered_mean, filtered_cov = args - # Get parameters and inputs for time index t - Q = _get_params(params.dynamics_covariance, 2, t) + # Get parameters and inputs for times t and t + 1 + Q = _get_params(params.dynamics_covariance, 2, t + 1) R = _get_params(params.emission_covariance, 2, t) - u = inputs[t] + u = inputs[t + 1] F_x = F(filtered_mean, u) # Prediction step @@ -272,9 +283,9 @@ def _step(carry, args): next_state = carry key, filtered_mean, filtered_cov, t = args - # Get parameters and inputs for time index t - Q = _get_params(params.dynamics_covariance, 2, t) - u = inputs[t] + # Get parameters and inputs for time index t + 1 + Q = _get_params(params.dynamics_covariance, 2, t + 1) + u = inputs[t + 1] # Condition on next state smoothed_mean, smoothed_cov = _condition_on(filtered_mean, filtered_cov, f, F, Q, u, next_state, 1) From 5c0e1206e5c6837a3399e3bf9a0cb3201defc7bb Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 14/24] test: check UKF input and parameter indexing, prediction lengths, and covariance lengths --- .../nonlinear_gaussian_ssm/inference_test.py | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/dynamax/nonlinear_gaussian_ssm/inference_test.py b/dynamax/nonlinear_gaussian_ssm/inference_test.py index 695a5ed23..744989c31 100644 --- a/dynamax/nonlinear_gaussian_ssm/inference_test.py +++ b/dynamax/nonlinear_gaussian_ssm/inference_test.py @@ -1,5 +1,5 @@ """ -Tests for the extended Kalman filter and smoother. +Tests for the extended and unscented Kalman filters and smoothers. """ import numpy as np import pytest @@ -10,6 +10,7 @@ from dynamax.linear_gaussian_ssm import lgssm_filter, lgssm_smoother, lgssm_posterior_sample from dynamax.linear_gaussian_ssm.inference_test_utils import assert_sample_moments from dynamax.nonlinear_gaussian_ssm.inference_ekf import extended_kalman_filter, extended_kalman_smoother, extended_kalman_posterior_sample +from dynamax.nonlinear_gaussian_ssm.inference_ukf import unscented_kalman_filter, unscented_kalman_smoother, UKFHyperParams from dynamax.nonlinear_gaussian_ssm.inference_test_utils import lgssm_to_nlgssm, random_lgssm_args, random_nlgssm_args, make_affine_timing_case from dynamax.nonlinear_gaussian_ssm.models import ParamsNLGSSM from dynamax.nonlinear_gaussian_ssm.sarkka_lib import ekf, eks @@ -138,12 +139,16 @@ def test_extended_kalman_sampler_nonlinear(key=0, num_timesteps=15, sample_size= "filter_fn,smoother_fn", [ (extended_kalman_filter, extended_kalman_smoother), + ( + partial(unscented_kalman_filter, hyperparams=UKFHyperParams()), + partial(unscented_kalman_smoother, hyperparams=UKFHyperParams()), + ), ], - ids=["extended"], + ids=["extended", "unscented"], ) @pytest.mark.parametrize("num_timesteps", [1, 3]) def test_kalman_inference_input_indexing(filter_fn, smoother_fn, num_timesteps): - """Check EKF input and parameter indexing by comparing filtering, smoothing, + """Check EKF and UKF input and parameter indexing by comparing filtering, smoothing, and prediction results with exact LGSSM results. """ params, emissions, inputs, expected = make_affine_timing_case(num_timesteps) @@ -158,13 +163,13 @@ def test_kalman_inference_input_indexing(filter_fn, smoother_fn, num_timesteps): @pytest.mark.parametrize( "filter_fn", - [extended_kalman_filter], - ids=["extended"], + [extended_kalman_filter, partial(unscented_kalman_filter, hyperparams=UKFHyperParams())], + ids=["extended", "unscented"], ) @pytest.mark.parametrize("num_timesteps", [1, 4]) @pytest.mark.parametrize("time_varying_covariance", [False, True]) def test_kalman_prediction_length_without_inputs(filter_fn, num_timesteps, time_varying_covariance): - """Check EKF prediction lengths and values without inputs, for static and + """Check EKF and UKF prediction lengths and values without inputs, for static and time-varying dynamics covariance. Check both full-output and prediction-only requests. @@ -207,12 +212,15 @@ def test_extended_kalman_posterior_sample_input_indexing(num_timesteps): def test_kalman_inference_rejects_short_dynamics_covariance(): - """Check that the EKF filter and smoother reject incorrect dynamics covariance lengths.""" + """Check that the EKF and UKF filters and the EKF smoother + reject incorrect dynamics covariance lengths. + """ params, emissions, inputs, _ = make_affine_timing_case() filtered = extended_kalman_filter(params, emissions, inputs=inputs) params = params._replace(dynamics_covariance=params.dynamics_covariance[:-1]) for inference in ( lambda y, u: extended_kalman_filter(params, y, inputs=u), + lambda y, u: unscented_kalman_filter(params, y, UKFHyperParams(), inputs=u), lambda y, u: extended_kalman_smoother(params, y, filtered_posterior=filtered, inputs=u), ): with pytest.raises(AssertionError): From f1aef25e3f7325c61a249b6f3b49314d98e09555 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 15/24] fix: correct UKF input and parameter indexing, prediction lengths, and state dimension calculation --- .../nonlinear_gaussian_ssm/inference_ukf.py | 30 ++++++++++++------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/dynamax/nonlinear_gaussian_ssm/inference_ukf.py b/dynamax/nonlinear_gaussian_ssm/inference_ukf.py index b67db3863..11371fa38 100644 --- a/dynamax/nonlinear_gaussian_ssm/inference_ukf.py +++ b/dynamax/nonlinear_gaussian_ssm/inference_ukf.py @@ -179,11 +179,15 @@ def unscented_kalman_filter(params: ParamsNLGSSM, inputs: optional array of inputs. Returns: - filtered_posterior: posterior object. + filtered_posterior: posterior object. Predicted entry ``t`` is for state ``t + 1``. + Predicted fields have length ``T`` when inputs are omitted and + dynamics covariance is static, otherwise ``T - 1``. """ num_timesteps = len(emissions) - state_dim = params.dynamics_covariance.shape[0] + predict_final = inputs is None and params.dynamics_covariance.ndim == 2 + assert params.dynamics_covariance.ndim != 3 or params.dynamics_covariance.shape[0] == num_timesteps + state_dim = params.initial_mean.shape[0] # Compute lambda and weights from from hyperparameters alpha, beta, kappa = hyperparams.alpha, hyperparams.beta, hyperparams.kappa @@ -195,12 +199,11 @@ def unscented_kalman_filter(params: ParamsNLGSSM, f, h = (_process_fn(fn, inputs) for fn in (f, h)) inputs = _process_input(inputs, num_timesteps) - def _step(carry, t): + def _step(carry, t, predict=True): """One step of the UKF""" ll, pred_mean, pred_cov = carry # Get parameters and inputs for time t - Q = _get_params(params.dynamics_covariance, 2, t) R = _get_params(params.emission_covariance, 2, t) u = inputs[t] y = emissions[t] @@ -214,7 +217,10 @@ def _step(carry, t): ll += log_likelihood # Predict the next state - pred_mean, pred_cov, _ = _predict(filtered_mean, filtered_cov, f, Q, lamb, w_mean, w_cov, u) + if predict: + Q = _get_params(params.dynamics_covariance, 2, t + 1) + u_next = inputs[t + 1] if not predict_final else u + pred_mean, pred_cov, _ = _predict(filtered_mean, filtered_cov, f, Q, lamb, w_mean, w_cov, u_next) # Build carry and output states carry = (ll, pred_mean, pred_cov) @@ -231,7 +237,11 @@ def _step(carry, t): # Run the Unscented Kalman Filter carry = (0.0, params.initial_mean, params.initial_covariance) - (ll, *_), outputs = lax.scan(_step, carry, jnp.arange(num_timesteps)) + carry, outputs = lax.scan(_step, carry, jnp.arange(num_timesteps - 1)) + (ll, _, _), final_outputs = _step(carry, num_timesteps - 1, predict=predict_final) + outputs = {key: jnp.concatenate((val, final_outputs[key][None])) for key, val in outputs.items()} + if not predict_final: + outputs = {key: val[:-1] if key.startswith("predicted") else val for key, val in outputs.items()} outputs = {"marginal_loglik": ll, **outputs} posterior_filtered = PosteriorGSSMFiltered( **outputs, @@ -257,7 +267,7 @@ def unscented_kalman_smoother(params: ParamsNLGSSM, """ num_timesteps = len(emissions) - state_dim = params.dynamics_covariance.shape[0] + state_dim = params.initial_mean.shape[0] # Run the unscented Kalman filter ukf_posterior = unscented_kalman_filter(params, emissions, hyperparams, inputs) @@ -281,10 +291,10 @@ def _step(carry, args): smoothed_mean_next, smoothed_cov_next = carry t, filtered_mean, filtered_cov = args - # Get parameters and inputs for time t - Q = _get_params(params.dynamics_covariance, 2, t) + # Get parameters and inputs for times t and t + 1 + Q = _get_params(params.dynamics_covariance, 2, t + 1) R = _get_params(params.emission_covariance, 2, t) - u = inputs[t] + u = inputs[t + 1] y = emissions[t] # Prediction step From a6c53a317df80a55cefb2b67f8fbb65554c9d8c4 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 16/24] test: check GGSSM input and parameter indexing and covariance lengths --- .../inference_test.py | 41 +++++++++++++++++-- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/dynamax/generalized_gaussian_ssm/inference_test.py b/dynamax/generalized_gaussian_ssm/inference_test.py index 47f2358aa..c357af809 100644 --- a/dynamax/generalized_gaussian_ssm/inference_test.py +++ b/dynamax/generalized_gaussian_ssm/inference_test.py @@ -1,13 +1,16 @@ """ Tests for inference in the generalized Gaussian SSM. """ +import numpy as np +import pytest import jax.numpy as jnp +from jax import jit from dynamax.generalized_gaussian_ssm.models import ParamsGGSSM -from dynamax.generalized_gaussian_ssm.inference import conditional_moments_gaussian_smoother, EKFIntegrals, UKFIntegrals +from dynamax.generalized_gaussian_ssm.inference import conditional_moments_gaussian_filter, conditional_moments_gaussian_smoother, EKFIntegrals, UKFIntegrals, GHKFIntegrals from dynamax.nonlinear_gaussian_ssm.inference_ekf import extended_kalman_smoother from dynamax.nonlinear_gaussian_ssm.inference_ukf import unscented_kalman_smoother, UKFHyperParams -from dynamax.nonlinear_gaussian_ssm.inference_test_utils import random_nlgssm_args +from dynamax.nonlinear_gaussian_ssm.inference_test_utils import random_nlgssm_args, make_affine_timing_case from dynamax.utils.utils import has_tpu from functools import partial @@ -71,4 +74,36 @@ def test_ukf(key=1, num_timesteps=15): assert allclose(ukf_post.filtered_means, ggf_post.filtered_means) assert allclose(ukf_post.filtered_covariances, ggf_post.filtered_covariances) assert allclose(ukf_post.smoothed_means, ggf_post.smoothed_means) - assert allclose(ukf_post.smoothed_covariances, ggf_post.smoothed_covariances) \ No newline at end of file + assert allclose(ukf_post.smoothed_covariances, ggf_post.smoothed_covariances) + +@pytest.mark.parametrize("num_timesteps", [1, 3]) +@pytest.mark.parametrize("integrals", [EKFIntegrals(), UKFIntegrals(), GHKFIntegrals(order=3)]) +def test_generalized_gaussian_inference_input_indexing(num_timesteps, integrals): + """Check GGSSM input and parameter indexing against exact LGSSM filtering and smoothing results. + + Also check that incorrect dynamics covariance lengths are rejected. + """ + params, emissions, inputs, expected = make_affine_timing_case(num_timesteps) + generalized_params = ParamsGGSSM( + initial_mean=params.initial_mean, + initial_covariance=params.initial_covariance, + dynamics_function=params.dynamics_function, + dynamics_covariance=params.dynamics_covariance, + emission_mean_function=params.emission_function, + emission_cov_function=lambda state, inpt: params.emission_covariance, + ) + posterior = jit(lambda y, u: conditional_moments_gaussian_smoother( + generalized_params, integrals, y, inputs=u + ))(emissions, inputs) + for field in posterior._fields: + value = getattr(posterior, field) + if value is not None: + np.testing.assert_allclose(value, expected[field], atol=2e-5, rtol=2e-5) + + # Check dynamics covariance lengths in both filtering and smoothing. + short_params = generalized_params._replace(dynamics_covariance=params.dynamics_covariance[:-1]) + with pytest.raises(AssertionError): + conditional_moments_gaussian_filter(short_params, integrals, emissions, inputs=inputs) + with pytest.raises(AssertionError): + conditional_moments_gaussian_smoother( + short_params, integrals, emissions, filtered_posterior=posterior, inputs=inputs) From ed53cb15a140dfea7ca3a844b0f11a2fc42d749a Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:36 -0700 Subject: [PATCH 17/24] fix: correct GGSSM input and parameter indexing --- dynamax/generalized_gaussian_ssm/inference.py | 33 +++++++++++-------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/dynamax/generalized_gaussian_ssm/inference.py b/dynamax/generalized_gaussian_ssm/inference.py index 324a4a111..8f1daf36e 100644 --- a/dynamax/generalized_gaussian_ssm/inference.py +++ b/dynamax/generalized_gaussian_ssm/inference.py @@ -183,11 +183,11 @@ def _predict(prior_mean: Float[Array, "state_dim"], p(x_{t+1}) = N(x_{t+1} | mu_pred, Sigma_pred) where mu_pred = gaussian_expec(f, m, P) - \approx \int f(x_t, u) N(x_t | m, P) dx_t - Sigma_pred = gaussian_crosscov(f(x_t, u_t), f(x_t, u_t); m, P) + Q - \approx \int (f(x_t, u) - mu_pred)(f(x_t, u) - mu_pred)^T N(x_t | m, P)dx_t + Q - cross_pred = gaussian_crosscov(x_t, f(x_t, u_t); m, P) - \approx \int (x_t - m)(f(x_t, u_t) - mu_pred)^T N(x_t | m, P)dx_t + \approx \int f(x_t, u_{t+1}) N(x_t | m, P) dx_t + Sigma_pred = gaussian_crosscov(f(x_t, u_{t+1}), f(x_t, u_{t+1}); m, P) + Q + \approx \int (f(x_t, u_{t+1}) - mu_pred)(f(x_t, u_{t+1}) - mu_pred)^T N(x_t | m, P)dx_t + Q + cross_pred = gaussian_crosscov(x_t, f(x_t, u_{t+1}); m, P) + \approx \int (x_t - m)(f(x_t, u_{t+1}) - mu_pred)^T N(x_t | m, P)dx_t Args: prior_mean (state_dim,): prior mean. @@ -232,7 +232,7 @@ def _condition_on(prior_mean: Float[Array, "state_dim"], Gaussian approximation. p(x_t | y_t, u_t, y_{1:t-1}, u_{1:t-1}) - propto p(x_t | y_{1:t-1}, u_{1:t-1}) p(y_t | x_t, u_t) + propto p(x_t | y_{1:t-1}, u_{1:t}) p(y_t | x_t, u_t) \approx N(x_t | m, P) ArbitraryDist(y_t | emission_mean(x_t, u_t), emission_cov(x_t, u_t)) \approx N(x_t | mu_cond, Sigma_cond) @@ -305,6 +305,7 @@ def conditional_moments_gaussian_filter(model_params: ParamsGGSSM, """ num_timesteps = len(emissions) + assert model_params.dynamics_covariance.ndim != 3 or model_params.dynamics_covariance.shape[0] == num_timesteps # Process dynamics function and conditional emission moments to take in control inputs f = model_params.dynamics_function @@ -319,12 +320,11 @@ def conditional_moments_gaussian_filter(model_params: ParamsGGSSM, # Emission distribution emission_dist = model_params.emission_dist - def _step(carry, t): + def _step(carry, t, predict=True): """One step of the CMGF""" ll, pred_mean, pred_cov = carry - # Get parameters and inputs for time index t - Q = _get_params(model_params.dynamics_covariance, 2, t) + # Get the input and observation at time t u = inputs[t] y = emissions[t] @@ -334,13 +334,17 @@ def _step(carry, t): ll += log_likelihood # Predict the next state - pred_mean, pred_cov, _ = _predict(filtered_mean, filtered_cov, f, Q, u, g_ev, g_cov) + if predict: + Q = _get_params(model_params.dynamics_covariance, 2, t + 1) + pred_mean, pred_cov, _ = _predict(filtered_mean, filtered_cov, f, Q, inputs[t + 1], g_ev, g_cov) return (ll, pred_mean, pred_cov), (filtered_mean, filtered_cov) # Run the general linearization filter carry = (0.0, model_params.initial_mean, model_params.initial_covariance) - (ll, _, _), (filtered_means, filtered_covs) = lax.scan(_step, carry, jnp.arange(num_timesteps)) + carry, outputs = lax.scan(_step, carry, jnp.arange(num_timesteps - 1)) + (ll, _, _), final_outputs = _step(carry, num_timesteps - 1, predict=False) + filtered_means, filtered_covs = [jnp.concatenate((a, b[None])) for a, b in zip(outputs, final_outputs)] return PosteriorGSSMFiltered(marginal_loglik=ll, filtered_means=filtered_means, filtered_covariances=filtered_covs) @@ -366,6 +370,7 @@ def conditional_moments_gaussian_smoother(model_params: ParamsGGSSM, """ num_timesteps = len(emissions) + assert model_params.dynamics_covariance.ndim != 3 or model_params.dynamics_covariance.shape[0] == num_timesteps # Get filtered posterior if filtered_posterior is None: @@ -386,9 +391,9 @@ def _step(carry, args): smoothed_mean_next, smoothed_cov_next = carry t, filtered_mean, filtered_cov = args - # Get parameters and inputs for time index t - Q = _get_params(model_params.dynamics_covariance, 2, t) - u = inputs[t] + # Get parameters and inputs for time index t + 1 + Q = _get_params(model_params.dynamics_covariance, 2, t + 1) + u = inputs[t + 1] # Prediction step pred_mean, pred_cov, pred_cross = _predict(filtered_mean, filtered_cov, f, Q, u, g_ev, g_cov) From 2538a1c8715c6e48be7f32070405d59ff6f94108 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Thu, 1 Oct 2026 16:29:37 -0700 Subject: [PATCH 18/24] docs: document the state indices and lengths of prediction arrays --- dynamax/linear_gaussian_ssm/inference.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/dynamax/linear_gaussian_ssm/inference.py b/dynamax/linear_gaussian_ssm/inference.py index 6e80f2ae0..6c98235ad 100644 --- a/dynamax/linear_gaussian_ssm/inference.py +++ b/dynamax/linear_gaussian_ssm/inference.py @@ -121,13 +121,15 @@ class PosteriorGSSMFiltered(NamedTuple): :param marginal_loglik: marginal log likelihood, $p(y_{1:T} \mid u_{1:T})$ :param filtered_means: array of filtered means $\mathbb{E}[z_t \mid y_{1:t}, u_{1:t}]$ :param filtered_covariances: array of filtered covariances $\mathrm{Cov}[z_t \mid y_{1:t}, u_{1:t}]$ + :param predicted_means: array of predicted means $\mathbb{E}[z_{t+1} \mid y_{1:t}, u_{1:t+1}]$ + :param predicted_covariances: array of predicted covariances $\mathrm{Cov}[z_{t+1} \mid y_{1:t}, u_{1:t+1}]$ """ marginal_loglik: Union[Scalar, Float[Array, " ntime"]] filtered_means: Optional[Float[Array, "ntime state_dim"]] = None filtered_covariances: Optional[Float[Array, "ntime state_dim state_dim"]] = None - predicted_means: Optional[Float[Array, "ntime state_dim"]] = None - predicted_covariances: Optional[Float[Array, "ntime state_dim state_dim"]] = None + predicted_means: Optional[Float[Array, "num_predictions state_dim"]] = None + predicted_covariances: Optional[Float[Array, "num_predictions state_dim state_dim"]] = None class PosteriorGSSMSmoothed(NamedTuple): From 977ebf271c41f8671ce2f521ac823fa1698996c9 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Fri, 2 Oct 2026 21:07:05 -0700 Subject: [PATCH 19/24] fix: correct transition indexing in HMM posterior sampling --- dynamax/hidden_markov_model/inference.py | 2 +- dynamax/hidden_markov_model/inference_test.py | 11 +++++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/dynamax/hidden_markov_model/inference.py b/dynamax/hidden_markov_model/inference.py index 5a8a95db2..21a3400b0 100644 --- a/dynamax/hidden_markov_model/inference.py +++ b/dynamax/hidden_markov_model/inference.py @@ -554,7 +554,7 @@ def _step(carry, args): keys = jr.split(key, num_timesteps) last_state = jr.choice(keys[-1], a=num_states, p=filtered_probs[-1]) _, states = lax.scan( - _step, last_state, (jnp.arange(1, num_timesteps), keys[:-1], filtered_probs[:-1]), + _step, last_state, (jnp.arange(num_timesteps - 1), keys[:-1], filtered_probs[:-1]), reverse=True ) diff --git a/dynamax/hidden_markov_model/inference_test.py b/dynamax/hidden_markov_model/inference_test.py index b8aaa2554..ca3fae2ba 100644 --- a/dynamax/hidden_markov_model/inference_test.py +++ b/dynamax/hidden_markov_model/inference_test.py @@ -366,3 +366,14 @@ def test_parallel_posterior_sample( # Compare the joint distributions assert jnp.allclose(blj_sample, blj, rtol=0, atol=eps) + + +def test_hmm_posterior_sample_nonstationary(): + """Check HMM posterior sampling transition indexing: + transitions[t] connects states[t] and states[t + 1]. + """ + transitions = jnp.stack([jnp.eye(2), jnp.fliplr(jnp.eye(2))]) + _, states = core.hmm_posterior_sample( + jr.PRNGKey(0), jnp.array([0.5, 0.5]), transitions, jnp.zeros((3, 2))) + assert states[0] == states[1] + assert states[1] != states[2] From 8a1c304818075a600f32714f7bed0bc963a9b3fb Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Fri, 2 Oct 2026 21:07:05 -0700 Subject: [PATCH 20/24] fix: correct transition indexing in HMM fixed-lag smoothing --- .../demos/fixed_lag_smoother_hmm.ipynb | 4 ++-- dynamax/hidden_markov_model/inference.py | 2 +- dynamax/hidden_markov_model/inference_test.py | 11 +++++++++++ 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb index 472b0023e..e4763a4d5 100644 --- a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb +++ b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb @@ -77,7 +77,7 @@ "\n", " # Get parameters for time t\n", " A_fwd = core._get_params(transition_matrix, 2, t_offset - 1)\n", - " A_bwd = core._get_params(transition_matrix, 2, t_offset)\n", + " A_bwd = core._get_params(transition_matrix, 2, t_offset - 1)\n", " ll = log_likelihoods[t_offset]\n", "\n", " # Shift window forward by 1\n", @@ -175,7 +175,7 @@ "\n", " # Get parameters for time t\n", " A_fwd = core._get_params(transition_matrix, 2, t - 1)\n", - " A_bwd = core._get_params(transition_matrix, 2, t)\n", + " A_bwd = core._get_params(transition_matrix, 2, t - 1)\n", " ll = log_likelihoods[t]\n", "\n", " # Shift window forward by 1\n", diff --git a/dynamax/hidden_markov_model/inference.py b/dynamax/hidden_markov_model/inference.py index 21a3400b0..47456f6c8 100644 --- a/dynamax/hidden_markov_model/inference.py +++ b/dynamax/hidden_markov_model/inference.py @@ -380,7 +380,7 @@ def _step(carry, t): # Get parameters for time t A_fwd = get_trans_mat(transition_matrix, transition_fn, t-1) - A_bwd = get_trans_mat(transition_matrix, transition_fn, t) + A_bwd = A_fwd ll = log_likelihoods[t] # Shift window forward by 1 diff --git a/dynamax/hidden_markov_model/inference_test.py b/dynamax/hidden_markov_model/inference_test.py index ca3fae2ba..c7c1ebbed 100644 --- a/dynamax/hidden_markov_model/inference_test.py +++ b/dynamax/hidden_markov_model/inference_test.py @@ -377,3 +377,14 @@ def test_hmm_posterior_sample_nonstationary(): jr.PRNGKey(0), jnp.array([0.5, 0.5]), transitions, jnp.zeros((3, 2))) assert states[0] == states[1] assert states[1] != states[2] + + +def test_hmm_fixed_lag_smoother_nonstationary(): + """Check HMM fixed-lag smoothing transition indexing: + both messages use transition index t - 1 at time t. + """ + transitions = jnp.stack([jnp.eye(2), jnp.fliplr(jnp.eye(2))]) + log_likelihoods = jnp.log(jnp.array([[1., 1.], [0.9, 0.1], [1., 1.]])) + posterior = core.hmm_fixed_lag_smoother( + jnp.array([0.5, 0.5]), transitions, log_likelihoods, window_size=2) + assert jnp.allclose(posterior.smoothed_probs[1], jnp.array([[0.9, 0.1], [0.9, 0.1]])) From 606a8b06e57d0838a546b5c19fd9fda1742e6c03 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Fri, 2 Oct 2026 21:07:05 -0700 Subject: [PATCH 21/24] fix: use the first input to compute HMM initial probabilities --- .../models/abstractions.py | 2 +- .../hidden_markov_model/models/test_models.py | 20 +++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/dynamax/hidden_markov_model/models/abstractions.py b/dynamax/hidden_markov_model/models/abstractions.py index 533b0ffcf..e819f2246 100644 --- a/dynamax/hidden_markov_model/models/abstractions.py +++ b/dynamax/hidden_markov_model/models/abstractions.py @@ -602,7 +602,7 @@ def _inference_args(self, params: HMMParameterSet, emissions: Array, inputs: Optional[Array]) -> Tuple: """Return the arguments needed for inference.""" - return (self.initial_component._compute_initial_probs(params.initial, inputs), + return (self.initial_component._compute_initial_probs(params.initial, pytree_slice(inputs, 0)), self.transition_component._compute_transition_matrices(params.transitions, inputs), self.emission_component._compute_conditional_logliks(params.emissions, emissions, inputs)) diff --git a/dynamax/hidden_markov_model/models/test_models.py b/dynamax/hidden_markov_model/models/test_models.py index ff1af284e..df174ba62 100644 --- a/dynamax/hidden_markov_model/models/test_models.py +++ b/dynamax/hidden_markov_model/models/test_models.py @@ -5,6 +5,9 @@ import jax.numpy as jnp import jax.random as jr import pytest +from dynamax.hidden_markov_model.models.abstractions import HMMInitialState +from dynamax.hidden_markov_model.models.initial import StandardHMMInitialState +from tensorflow_probability.substrates.jax import distributions as tfd from jax import vmap from jax.tree_util import tree_leaves, tree_map @@ -369,3 +372,20 @@ def init_kmeans(key, emissions, inputs): vmapped_params = vmap(init_kmeans)(batch_keys, batch_emissions, batch_inputs) assert vmapped_params.emissions.biases.shape == (batch_size, num_states) assert jnp.all(jnp.isfinite(vmapped_params.emissions.biases)) + + +def test_hmm_input_dependent_initial_distribution(): + """Check that HMM inference uses the first input to compute initial probabilities.""" + class InputInitialState(StandardHMMInitialState): + _compute_initial_probs = HMMInitialState._compute_initial_probs + + def distribution(self, params, inputs=None): + return tfd.Categorical(logits=inputs) + + hmm = models.GaussianHMM(num_states=2, emission_dim=1) + params, _ = hmm.initialize(jr.PRNGKey(0)) + hmm.initial_component = InputInitialState(2) + inputs = jnp.array([[3., 0.], [0., 3.], [2., 1.]]) + posterior = hmm.filter(params, jnp.zeros((3, 1)), inputs=inputs) + expected = hmm.initial_distribution(params, inputs[0]).probs_parameter() + assert jnp.allclose(posterior.predicted_probs[0], expected) From c66c3d46f7c122966e7a65249e9fdf36a788e40f Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Sat, 3 Oct 2026 12:10:20 -0700 Subject: [PATCH 22/24] fix: add missing initial_probs in the HMM fixed-lag notebook --- .../hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb index e4763a4d5..fd21e3e67 100644 --- a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb +++ b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb @@ -114,6 +114,7 @@ " filtered_probs=filtered_probs,\n", " predicted_probs=predicted_probs,\n", " smoothed_probs=smoothed_probs,\n", + " initial_probs=smoothed_probs[0],\n", " )\n", "\n", " return (log_normalizers, filtered_probs, predicted_probs), post\n", @@ -151,6 +152,7 @@ " filtered_probs=filtered_probs,\n", " predicted_probs=predicted_probs,\n", " smoothed_probs=smoothed_probs,\n", + " initial_probs=smoothed_probs[0],\n", " )\n", "\n", " return posts" @@ -217,6 +219,7 @@ " filtered_probs=filtered_probs,\n", " predicted_probs=predicted_probs,\n", " smoothed_probs=smoothed_probs,\n", + " initial_probs=smoothed_probs[0],\n", " )\n", "\n", " return (log_normalizers, filtered_probs, predicted_probs, bmatrices), post\n", @@ -245,6 +248,7 @@ " filtered_probs=filtered_probs,\n", " predicted_probs=predicted_probs,\n", " smoothed_probs=smoothed_probs,\n", + " initial_probs=smoothed_probs[0],\n", " )\n", "\n", " return posts" From 4bd43c67a5965dfcbf27ff086249c0ce5cd70bd4 Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Sat, 3 Oct 2026 12:30:35 -0700 Subject: [PATCH 23/24] fix: wait for computations before recording HMM fixed-lag benchmark times --- .../hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb index fd21e3e67..d7e863f0c 100644 --- a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb +++ b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb @@ -45,7 +45,7 @@ "import jax.numpy as jnp\n", "import jax.random as jr\n", "import jax.lax as lax\n", - "from jax import vmap" + "from jax import vmap, block_until_ready" ] }, { @@ -423,12 +423,12 @@ " for window_size in window_grid:\n", " print(f\"Window of size: {window_size}\")\n", " start = time.time()\n", - " _ = hmm_fixed_lag_smoother_iterative(*args, window_size)\n", + " _ = block_until_ready(hmm_fixed_lag_smoother_iterative(*args, window_size))\n", " it_time = time.time() - start\n", " it_times.append(it_time)\n", " print(f\"Iterative version took {it_time} seconds.\")\n", " start = time.time()\n", - " _ = hmm_fixed_lag_smoother_vectorized(*args, window_size)\n", + " _ = block_until_ready(hmm_fixed_lag_smoother_vectorized(*args, window_size))\n", " vec_time = time.time() - start\n", " vec_times.append(vec_time)\n", " print(f\"Vectorized version took {vec_time} seconds.\")\n", From cbfd6d7f8692236d923175766935bef22d2841af Mon Sep 17 00:00:00 2001 From: nkbranigan <81596699+nkbranigan@users.noreply.github.com> Date: Sat, 3 Oct 2026 12:30:35 -0700 Subject: [PATCH 24/24] docs: update HMM fixed-lag notebook outputs after rerunning all cells --- .../demos/fixed_lag_smoother_hmm.ipynb | 86 +++++++++---------- 1 file changed, 42 insertions(+), 44 deletions(-) diff --git a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb index d7e863f0c..0467d9e2b 100644 --- a/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb +++ b/dynamax/hidden_markov_model/demos/fixed_lag_smoother_hmm.ipynb @@ -336,74 +336,72 @@ "output_type": "stream", "text": [ "Window of size: 2\n", - "Iterative version took 0.749182939529419 seconds.\n", - "Vectorized version took 0.3245062828063965 seconds.\n", + "Iterative version took 0.36290621757507324 seconds.\n", + "Vectorized version took 0.15166306495666504 seconds.\n", "Window of size: 500\n", - "Iterative version took 0.9897620677947998 seconds.\n", - "Vectorized version took 0.5744988918304443 seconds.\n", + "Iterative version took 0.383882999420166 seconds.\n", + "Vectorized version took 0.24146676063537598 seconds.\n", "Window of size: 1000\n", - "Iterative version took 1.2942566871643066 seconds.\n", - "Vectorized version took 0.7614669799804688 seconds.\n", + "Iterative version took 0.6629061698913574 seconds.\n", + "Vectorized version took 0.2890188694000244 seconds.\n", "Window of size: 1500\n", - "Iterative version took 1.7407751083374023 seconds.\n", - "Vectorized version took 1.0755679607391357 seconds.\n", + "Iterative version took 0.9182291030883789 seconds.\n", + "Vectorized version took 0.34940099716186523 seconds.\n", "Window of size: 2000\n", - "Iterative version took 2.0848448276519775 seconds.\n", - "Vectorized version took 1.2289578914642334 seconds.\n", + "Iterative version took 1.1316211223602295 seconds.\n", + "Vectorized version took 0.45354294776916504 seconds.\n", "Window of size: 2500\n", - "Iterative version took 2.5447559356689453 seconds.\n", - "Vectorized version took 1.5629408359527588 seconds.\n", + "Iterative version took 1.3302099704742432 seconds.\n", + "Vectorized version took 0.5359401702880859 seconds.\n", "Window of size: 3000\n", - "Iterative version took 3.1088812351226807 seconds.\n", - "Vectorized version took 1.8927171230316162 seconds.\n", + "Iterative version took 1.6161599159240723 seconds.\n", + "Vectorized version took 0.6290450096130371 seconds.\n", "Window of size: 3500\n", - "Iterative version took 3.738128185272217 seconds.\n", - "Vectorized version took 2.3711469173431396 seconds.\n", + "Iterative version took 1.788140058517456 seconds.\n", + "Vectorized version took 0.7084269523620605 seconds.\n", "Window of size: 4000\n", - "Iterative version took 4.1995158195495605 seconds.\n", - "Vectorized version took 2.3815088272094727 seconds.\n", + "Iterative version took 2.0038373470306396 seconds.\n", + "Vectorized version took 0.7731838226318359 seconds.\n", "Window of size: 4500\n", - "Iterative version took 5.630947828292847 seconds.\n", - "Vectorized version took 3.7188751697540283 seconds.\n", + "Iterative version took 2.307873010635376 seconds.\n", + "Vectorized version took 1.0258479118347168 seconds.\n", "Window of size: 5000\n", - "Iterative version took 6.471923112869263 seconds.\n", - "Vectorized version took 4.330507040023804 seconds.\n", + "Iterative version took 2.236435890197754 seconds.\n", + "Vectorized version took 0.9173357486724854 seconds.\n", "Window of size: 5500\n", - "Iterative version took 5.676476955413818 seconds.\n", - "Vectorized version took 4.0590901374816895 seconds.\n", + "Iterative version took 2.602116823196411 seconds.\n", + "Vectorized version took 1.071423053741455 seconds.\n", "Window of size: 6000\n", - "Iterative version took 6.0062477588653564 seconds.\n", - "Vectorized version took 5.020081996917725 seconds.\n", + "Iterative version took 2.730694055557251 seconds.\n", + "Vectorized version took 1.270918369293213 seconds.\n", "Window of size: 6500\n", - "Iterative version took 6.005163908004761 seconds.\n", - "Vectorized version took 4.315921068191528 seconds.\n", + "Iterative version took 2.9723217487335205 seconds.\n", + "Vectorized version took 1.4774446487426758 seconds.\n", "Window of size: 7000\n", - "Iterative version took 6.611321926116943 seconds.\n", - "Vectorized version took 4.7442381381988525 seconds.\n", + "Iterative version took 3.1633460521698 seconds.\n", + "Vectorized version took 1.7359628677368164 seconds.\n", "Window of size: 7500\n", - "Iterative version took 7.334857940673828 seconds.\n", - "Vectorized version took 5.0066869258880615 seconds.\n", + "Iterative version took 3.680619955062866 seconds.\n", + "Vectorized version took 2.0683720111846924 seconds.\n", "Window of size: 8000\n", - "Iterative version took 7.030580043792725 seconds.\n", - "Vectorized version took 5.242753028869629 seconds.\n", + "Iterative version took 4.33176326751709 seconds.\n", + "Vectorized version took 2.1112937927246094 seconds.\n", "Window of size: 8500\n", - "Iterative version took 7.414398908615112 seconds.\n", - "Vectorized version took 5.604051113128662 seconds.\n", + "Iterative version took 4.318212985992432 seconds.\n", + "Vectorized version took 1.893186092376709 seconds.\n", "Window of size: 9000\n", - "Iterative version took 7.887447118759155 seconds.\n", - "Vectorized version took 5.811422109603882 seconds.\n" + "Iterative version took 4.430753946304321 seconds.\n", + "Vectorized version took 2.271007776260376 seconds.\n" ] }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+GElEQVR4nO3dd3hUVfrA8e8hJIQQSiChhhB6CTWEXhUEKYIKSBMFRQQbrj/d1V0V2dVdO4ooCAKKoIAoRVERlBI6CR0ChBZIaCEhpJCe8/vjDhiQMoS5cyeT9/M8eZi5M3PPO5fJOyfnnvsepbVGCCGE+ylmdQBCCCHMIQleCCHclCR4IYRwU5LghRDCTUmCF0IIN1Xc6gDy8/f318HBwVaHIYQQhUZkZOR5rXXA9R5zqQQfHBxMRESE1WEIIUShoZSKudFjMkQjhBBuShK8EEK4KVMTvFLqb0qpfUqpvUqpb5VS3ma2J4QQ4k+mjcErpaoBzwGNtNbpSqmFwBDgy9vZT3Z2NrGxsWRkZJgQZdHm7e1NYGAgnp6eVocihDCB2SdZiwMllVLZgA9w6nZ3EBsbS+nSpQkODkYp5fAAiyqtNQkJCcTGxlKzZk2rwxFCmMC0IRqtdRzwPnACOA1c1Fr/du3zlFJjlFIRSqmI+Pj4v+wnIyODChUqSHJ3MKUUFSpUkL+MhHBjpiV4pZQf0B+oCVQFSimlHr72eVrr6VrrMK11WEDAdadySnI3iRxXIdybmSdZuwPHtNbxWuts4AegvYntCSFEoRMZk8jna4+Ysm8zE/wJoK1SykcZXcVuQJSJ7ZmmfXvje+n48eN88803Dt33f//73+u2JYRwb7l5mk9XH+ahzzfzzdYTpGXmOLwNM8fgtwCLgO3AHltb081qz0wbN24ECpbgc3Ju/p92bYK/3JYQwn2dTc5gxMwtvLfiIL2bVOHHZztSqoTj57yYOg9eaz1Ba91Aa91Yaz1Ca51pZntm8fX1BeDll18mPDyc5s2bM2nSJHJzc3nppZdo1aoVTZs25fPPPwdgzZo1dOrUiX79+tGoUSMA7r//flq2bElISAjTp0+/sr/09HSaN2/O8OHDr2pryJAhLF++/EoMI0eOZNGiRTdsUwhROPxx4Cy9Pg5nx4kk3h3QlMlDmlPG25ypyi5Vi+ZWJv64j/2nkh26z0ZVyzDhvhC7nvv222/z/vvv89NPPwEwffp0ypYty7Zt28jMzKRDhw706NEDgO3bt7N3794rUxBnzZpF+fLlSU9Pp1WrVgwYMIC3336bKVOmsHPnzr+0NXjwYBYuXEifPn3Iysri999/Z+rUqcycOfO6bcpURyFcW2ZOLu/+epCZ64/RoHJppgwLpU5FX1PbLFQJ3tX89ttv7N69m0WLFgFw8eJFoqOj8fLyonXr1lcl3cmTJ7N48WIATp48SXR0NBUqVLjhvnv16sX48ePJzMzk119/pXPnzpQsWfKGbUqCF8J1HTufxrPfbmdvXDKPtqvBK70b4u3pYXq7hSrB29vTdhatNZ988gk9e/a8avuaNWsoVarUVfdXrVrFpk2b8PHxoWvXrrecf+7t7U3Xrl1ZsWIFCxYsYMiQITdtUwjhmhbviOXVxXvxLF6M6SNa0iOkstPalmJjt6F06dKkpKRcud+zZ0+mTp1KdnY2AIcOHSItLe0vr7t48SJ+fn74+Phw4MABNm/efOUxT0/PK6+/1uDBg5k9ezbh4eHce++9t9WmEMJaaZk5vLBwJ39bsIuQqmX5+blOTk3uUMh68FZr2rQpHh4eNGvWjJEjRzJ+/HiOHz9OaGgoWmsCAgJYsmTJX1537733Mm3aNBo2bEj9+vVp27btlcfGjBlD06ZNCQ0NZd68eVe9rkePHowYMYL+/fvj5eUFwOjRo+1qUwhhnb1xF3n22x3EJKQxvltdnr27DsU9nN+fVlprpzd6I2FhYfraBT+ioqJo2LChRRG5Pzm+QjiO1ppZG47z9i9RVChVgo+GNKdtrRufa3MEpVSk1jrseo9JD14IIRwgITWTlxbt5o8D5+jesCLvDWyGXykvS2OSBC+EEHdo45HzPD9/J0mXsnnjvkY82t41qt9KghdCiALKyc3j49+jmbL6MDX9SzF7VCtCqpa1OqwrJMELIUQBbD2WyOtL93LgTAoDWwYysV+IKeUG7oRrRSOEEHbKzdN4FHP+MMjZ5Az+93MUS3aeompZb6Y9HMq9jas4PQ57SIIXQhQ6czYd582foujTtApju9SmfuXSpreZlZPHlxuP8fGqaLJzNc/eXYdxXWvj4+W6adR1I3MjS5YsoV69elcKj9lr2rRp+Pj48Mgjj9xR+8HBwURERODv739H+xHCFfy85zQTlu2jYeUyrNh3hsU74ujesCLjutahZQ0/U9oMj47njWX7OBKfRrcGFXmtbyOC/Uvd+oUWkwTvBEuWLKFv3763leBzcnIYO3asiVEJUfhsOZrA8wt20jLIj7mj25CelcucTTF8ufEYA6ZupHXN8ozrWpuu9QIcMosl9sIl3loexS97z1Cjgg+zRoZxd4NKDngnziEJ/hZefvllqlevztNPPw3AG2+8ga+vL1prFi5cSGZmJg888AATJ04EYM6cObz//vsopWjatCnjxo1j2bJlrF27ljfffJPvv/+elJQUxo4dy6VLl6hduzazZs3Cz8+Prl270rx5c9avX8/QoUNJSUnB19eXYcOG0bt37ysx7dmzh6NHj+Lj48PYsWM5ceIEAB999BEdOnQgISGBoUOHEhcXR7t27XCli9mEKKiDZ1IYPSeCoPI+fPFoGN6eHnh7ejC+e12e6FyT+VtPMiP8KKNmb6NhlTKM61qb3o0rF+gK0ozsXKavO8pnaw4D8FLP+jzesaZTCoQ5UuFK8L+8DGf2OHaflZtAr7dv+PDgwYN5/vnnryT4hQsX8o9//IMNGzawdetWtNb069ePdevWUaFCBd588002btyIv78/iYmJlC9fnn79+tG3b18GDhwIGCUPPvnkE7p06cLrr7/OxIkT+eijjwDIysri8tW8b7zxBgBVq1a9UlL4008/Ze3atdSoUYNhw4bxt7/9jY4dO3LixAl69uxJVFQUEydOpGPHjrz++ussX76cmTNnOvaYCeFkp5LSeXTWVny8PPjqsdaU87n6AiIfr+I81rEmD7etwbJdp5i29gjPfbuDDyr4MKZzLQaEBtqVnLXW/B51jn//tJ8TiZfo06QK/+zTkGrlSpr11kxVuBK8BVq0aMG5c+c4deoU8fHx+Pn5sWfPHn777TdatGgBQGpqKtHR0ezatYtBgwZdGesuX778X/Z38eJFkpKS6NKlCwCPPvoogwYNuvL44MGDbxjLhg0bmDFjBuvXrwdg1apV7N+//8rjycnJpKamsm7dOn744QcA+vTpg5+fOeOSQjjDxUvZPDprK2mZOXw3rt1Nk61X8WIMbBnIgy2qsTLqLJ+tOcK/Fu/lo1XRPN6xJsPbBFH6BotrHDufxr9/3Mfqg/HUqejLvNFt6FCncJ+3Mi3BK6XqAwvybaoFvK61/qjAO71JT9tMgwYNYtGiRZw5c4bBgwcTExPDK6+8wpNPPnnV8z755JM7bit/meH8Tp8+zeOPP86yZcuurPqUl5fH5s2b8fb2vuN2hXBFGdm5PDEngpiES3z5WCsaVC5j1+uKFVP0DKlMj0aV2HQ0galrjvD2Lwf4dPVhHmlXg1EdauLvWwKAS1k5TPnjMF+EH8OreDFe7dOQR9sH42lBcTBHM3NN1oNa6+Za6+ZAS+ASsNis9sw0ePBg5s+fz6JFixg0aBA9e/Zk1qxZpKamAhAXF8e5c+e4++67+e6770hISAAgMTERuLrMcNmyZfHz8yM8PByAr7/++kpv/kays7MZNGgQ77zzDvXq1buyvUePHld9qVwexuncufOVtWN/+eUXLly44ICjIIRz5eZpnp+/k63HE/ngoWa0r337vWmlFO1r+/P142348ZmOdK4bwGdrjtDh7T94bcleFkacpNsHa/lszRH6NqvCH//XhdGdarlFcgfnDdF0A45orWOc1J5DhYSEkJKSQrVq1ahSpQpVqlQhKiqKdu3aAcY6qnPnziUkJIR//etfdOnSBQ8PD1q0aMGXX37JkCFDeOKJJ5g8eTKLFi3iq6++unKStVatWsyePfum7W/cuJGIiAgmTJjAhAkTAPj555+ZPHkyTz/9NE2bNiUnJ4fOnTszbdo0JkyYwNChQwkJCaF9+/YEBQWZfoyEcCStNRN/3Mev+87wWt9G3Nes6h3vs0lgWT4dHsrR+FSmrzvK/G0n+HqzplGVMnwytAVhwX8dUi3snFIuWCk1C9iutZ5yncfGAGMAgoKCWsbEXP0dIOVszSXHV7iiT1cf5r0VBxnTuRb/7G3O5/PMxQwOnU2hQx1/S66IdZSblQs2/e8QpZQX0A/47nqPa62na63DtNZhAQEBZocjhHBxiyJjeW/FQe5vXpWX721gWjuVy3rTuV5AoU7ut+KMgaZeGL33s05oSwhRiK05eI5/fL+bjnX8eXdgM4q5cfJ1Bmck+KHAt3eyA7lQxxxyXIUr2XUyiafmbad+pdJMfTgUr+LucaLTSqYeQaVUKeAe4IeC7sPb25uEhARJRg6mtSYhIUGmWAqXcPx8Go99uY3ypbz48rFWN5yrLm6PqbNotNZpwB0tSBgYGEhsbCzx8fEOikpc5u3tTWBgoNVhiCLufGomj87eSp7WzHmsNRVLS6fDUVz+SlZPT09q1qxpdRhCCBOkZebw2JfbOJucwbdPtKVWgK/VIbkVl0/wQgj3lJ2bx7h529l3KpnpI1rSIkhKajianMUQQjid1pp/fL+bdYfieev+xnRrWHhK8BYmkuCFEE733oqD/LA9jr91r8eQ1nKltVlkiEYI4TRaayb/fpjP1hxhaOsgnutWx+qQ3JokeCGEU2iteXN5FDPXH2NAaCD/6R/ikFWXxI1JghdCmC43T/PKD7tZGBHLyPbBvN63kVyl6gQyBi+EE6zcf5ZDZ1OsDsMSmTm5PPvtdhZGxDK+W10m3CfJ3VmkBy+EyVbuP8sTcyKoUMqLpc90INDPx+qQnOZSVg5Pfh1JePR5Xu3TkNGdalkdUpEiPXghTHQi4RIvLNxJ/UqlycrJY8ycSNKzcq0OyykupmczYuZWNhw+z7sDmkpyt4AkeCFMkpGdy7h5kSjgi0fDmDy0BVFnknlp0S63r610PjWTodM3szs2iU+HhfJQq+pWh1QkSYIXwiRvLNvHvlPJTBrcnOrlfbirQUX+3rMBP+0+zdS1R6wOzzRxSek8NG0TR8+n8sWjrejVpIrVIRVZMgYvhAm+izjJ/G0neapr7auu0hzbpRZRp5N5b8VB6lcq7XZXcB6JT2XEF1tIycxh7uNt3HIZvMJEevBCONj+U8m8umQv7WpV4IV76l31mFKKdwY0JaRqGcbP38nhc+4zs2bfqYs8NG0TWbl5zB/TVpK7C5AEL4QDJWdk89S8SMqW9GTy0BYU9/jrr1hJLw8+HxGGt2cxnpgTycX0bAsidayI44kMmb6ZEsWLsfDJdoRULWt1SAJJ8EI4jNaav3+3m5MX0vl0eCgBpUvc8LnVypVk6sMtib1wiee+3UFuXuE96br2UDwPz9xCgG8JvhvXXkr+uhBJ8EI4yBfhx/h13xle6dWAVnYMT7QKLs+/+zdm7aF43l1xwAkROt7Pe04z+qtt1PL3ZeHYdlQrV9LqkEQ+Zi/ZV04ptUgpdUApFaWUamdme0JYZeuxRN7+9QC9Glfm8Y72L1AztHUQI9rW4PO1R1myI87ECB1v4baTPPPNdpoFluPbMW3x973xXyzCGmbPovkY+FVrPVAp5QUUnUv4RJFxLiWDZ77ZTlB5H94d2PS2C2i9fl8jDp5N4R/f76ZWQCmaBpYzJ1AH+iL8KG8uj6JTXX8+H9ESHy+ZkOeKlFkXXCilygI7gVrazkbCwsJ0RESEKfEIYYac3DwenrmFnSeTWPJ0BxpULlOg/SSkZtJvygbytGbpMx1MWZc0L0+z8UgCSelZV7YpjC+jy99J+b+a/vyeuvo5kTEXmL7uKL0aV+ajIc0pUdzD4bEK+ymlIrXWYdd7zMyv3ZpAPDBbKdUMiATG2xbiFsItfLDyEJuPJvLBoGYFTu4AFXxLMP2Rlgycuolxc7fzzRNtHJY4tdas2HeGSSujOeiggmeDWgbyvwebXHeWkHAdZvbgw4DNQAet9Ral1MdAstb6tWueNwYYAxAUFNQyJibGlHiEcLTLRcSGtg7ifw82ccg+l+8+zdPfbGdIq+r878Emd1QvXWvN71HnmLTqEPtOJVMroBTju9WlURXji0hfeV6+19i2Xt52vcdKFC9G7QBfqeXuIqzqwccCsVrrLbb7i4CXr32S1no6MB2MIRoT4xHCYS4XEWtcrQwT7mvksP32aVqFqNN1mLL6MI2qluGRdsG3vQ+tNWsPxTNp5SF2xV6kRgUfPnyoGf2aVZUedxFjWoLXWp9RSp1UStXXWh8EugH7zWpPCGe5XESsmFJMHd4Sb0/HjkG/cE89DpxJZuKP+6lbsTTtalew63VaG2PsH648RGTMBQL9SvLugKY8EFoNT0nsRZLZp76fBebZZtAcBUaZ3J4QprtcRGzWyDCql3f8xLBixRSTBjfngc828tS8SJY90/GW7Ww5msAHKw+x9VgiVcp689YDjRnUsjpexSWxu7zUeDh/EII7OnzXpiZ4rfVO4LpjQ0IURpeLiD19V23ubmBeobDS3p7MeCSM/lPW88ScCH54qv11pyJGxiTy4cpDbDicQMXSJZjYL4QhravLzBZXl3QSDvwEUT/CiU1Qogy8dBg8PB3ajExeFcJOl4uIta9dgRfuqW96ezX9S/HJsFBGzd7Ki9/t4tNhoVdObO48mcSklYdYeygef18vXuvbiOFtghw+XCQc6Hw0RC0zkvqpHca2iiHQ+e/Q8D4o5vh0LAleCDtcLiJWzscoIubhpDVFu9QL4JVeDXnr5yim/HGYuxpU5KNVh1gVdQ4/H09e6dWAEe1qyIVGrkhrOLPbSOhRP0K8rRxFtTDoPtFI6hVqmxqCfCqEuAWtNS8u3MXJC+kssOCS/NGdarL/dDIfrDzEBysPUbakJy/1rM+j7YPxLSG/wi4lLw9it9qS+jJIOgGqGNToAGGPQ4M+ULaa08KRT4cQtzBz/TF+23+WV/s0tKTGuVLKmBMPBFXw4bGONSnj7dixWnEHcrPh+HojqR/4CVLPgocX1LrLGH6p3wtK+VsSmiR4IW7iYno2k1YeonvDirdVRMzRvD09+HBwc8vaF9eRkQwbPoJtMyEjCTx9oO490LCf8a+39TXxJcELcRPfbj1BWlYuf7unnly5KQy5ObD9K1jzP0iLh0b9oelgqH03eLpWuWRJ8ELcQFZOHrM3HKNDnQqyQpEwTppGr4SVrxknTIPaw7AFUK2l1ZHdkCR4IW7gx12nOJucyTsDmlodirDamT3w26twdA2UrwWD5xknTF38rzpJ8EJch9aaGeFHqV+pNF3qBVgdjrBK8mlY/SbsmAcly8G970DYY1Dcy+rI7CIJXojrWH/4PAfOpBRoAQ/hBrLSYOMnsOFjY5ZMu6eh84tQ0s/qyG6LJHghrmP6uqMElC5B/+ZVrQ5FOFNeLuz8Bv54E1LPQKP7ofsbUN66GVR3QhK8ENeIOp1MePR5XupZX2q6FCVHVsNvr8HZPRDYCh6aA0FtrI7qjtwywSulSgHpWus8pVQ9oAHwi9Y62/TohLDAF+HH8PHyYHibIKtDEc5w7oAxMyb6NygXBANnQciDLn8C1R729ODXAZ2UUn7Ab8A2YDAw3MzAhLDCmYsZLNsVx/A2NSjnUzhOpIkCSjwGGydD5Ffg5Qv3/BtaPwmejl8P1yr2JHiltb6klHoc+Exr/a5SaqfJcQlhiS83Hic3T/NYh8I55ipuIf0C7FsMuxbAyc1GBcdWj0OXl6GUfQurFCZ2JXilVDuMHvvjtm0yMCncTmpmDt9siaFX4yoEVXD8Qh7CIjlZcHgV7PoWDv0KuVngXx+6vW5cgVo20OoITWNPgh8PvAIs1lrvU0rVAlabG5YQzrdw20mSM3IY3Ul674We1hC33Ujqe7+H9ETw8TfmsDcbAlWau8UY+63cMsFrrddhjMNfvn8UeM7MoIRwtpzcPGauP0arYD9aBBWuuc4in6QTsHsB7JoPCYfBowQ06A3Nhhq1Yhy8YpKrM3WapFLqOJAC5AI5WmtZvk+4pF/2niEuKZ0J9zWyOhRxuzIuwv6lxrh6zHpjW40O0GG8UQjMBao6WsUZ8+Dv0lqfd0I7QhSI1povwo9S078U3Ruat86q5fLyIGIm1OleaC/cuSI3B478YQzBHPwZcjKgQh24+1Vo8hD41bA6QpcgFzqJIm/rsUR2xV7kzfsbU8xJS/FZYss0WPEKlK0Oj61w6spCDpN8CrbPMaY2ppyCkuWhxQhjCKZaaJEYV78d9lzoFAA8AQTnf77W+jE79q+B35RSGvhcaz39OvsfA4wBCAqSC0uE880IP0r5Ul4MCHXf2RScPwy/T4TqbeHsPpj7IIz6BXycv0LVbcvLg6N/QMRsOPgL6Fyo3Q16vwt1exaawl9WsKcHvxQIB1ZhjKXfjo5a6zilVEVgpVLqgO2k7RW2pD8dICwsTN/m/oW4I0fiU1kVdY7x3epS0stNZ//m5cKScVDcGx76Cs5Hw9wBMG8QPLIUSvhaHeH1pZ6DHXMh8ktIijFmwXR4DkIfLfxDTE5iT4L30Vr/oyA711rH2f49p5RaDLQm34wcIaz2RfgxShQvxoh2bjxmu2mKsRD0g19A6crGz8BZsHCE8TN0gev0grU21jeNmGWscZqXDcGdoPsEaNAXijt3wfPCzp4E/5NSqrfW+ufb2bGthk0xrXWK7XYP4N8FCVIIM5xPzeT77bEMCA3E39dNE8e5A/DHW9DwPmgy8M/tDfvCfZNh2TOwZKyR/IsVsy7OS4nGCdOI2ZAQbcx8af0EtBwFAfWsi6uQs/dCp38qpTKBbEABWmtd5havqwQsttXSLg58o7X+9U6CFcKR5myKISsnz30vbMrNMYZmSvhCn0l/PQEZOgIuJcCqCcbJyt7vOfckpdYQu83ore9bbMyECWwN90+FkAdcbn3TwsieC51KF2THtguimhXktUKYLT0rl7mbY+jesBK1A1x0DPpObfgITm2HQV+C7w1Wper4PFw6byxuUcofur5sflw5WbDjayOxn90LXqWh+XAIGwWVm5jffhFywwSvlGqgtT6glAq93uNa6+3mhSWEub7fHktiWhZjOteyOhRznN0Ha942yt6GPHDz597zH2OIZM3/jJ58mzHmxXVqJyx92kjslZvCfR9D44Gue6K3kLtZD/4FjOmLH1znMQ3cbUpEQpgsN08zc/0xmgWWpVWwG5YlyM2GxWONNUR7v3/r5ytljMenX4Bf/m5Mncw/Xu8IOZmw9l1YPwlKBcCQb6F+L5m3brIbJnit9Rjbv3c5LxwhzLcq6izHzqcxZVgL91xvNfwDOLMbBs+zvwSuR3FjZs3cAbD4SePLoU53x8QTFwlLnob4KGMopudbhW5t08LKwtPmQljji/CjBPqV5N6QylaH4nindsK694wyuA373t5rPUvC0G+hYkNYMAJObruzWLIzYOUE+KK7US9m2Hdw/2eS3J1IErwoUrafuMC24xd4rENNinu42cc/JxOWPGVcEHTv2wXbh3dZePgH8K0E8wbCuaiC7efkNvi8k3Git/lweHoz1OtRsH2JAnOzT7gQN/dF+FHKeBfnoVbVrQ7F8da+C+f2Qb/Jd1aCwLcijFhsXFT09YNGCV57ZafDin/BrB6QdQke/h76TynSFR2tdMsErwwPK6Vet90PUkq1Nj80IRzrRMIlft17huFta+Bbws3q7MVFGicwmz8M9Xre+f7K1zR68tlp8PUDkGZHQdgTm2FaR+PK2dBH4KlNjhvHFwViTw/+M6AdMNR2PwX41LSIhDDJrA3H8CimGNk+2OpQHCs7AxaPM0oQ9HzLcfut3BiGLYSLccbJ14zk6z8vKw1+eRlm3WvMcX9kqTH90ftW10IKs9mT4NtorZ8GMgC01hcAFylcIYR9ki5lsWDbSfo1q0alMt5Wh+NYq9+C8weNoZmS5Ry776C2RoGyM3tg/jDjyyS/4xtgagfYMtVYvPqpjVCrq2NjEAVmT4LPVkp5YMx9v1w+OM/UqIRwsHlbTpCencsTnd2sLMGJLcZVqC1HmjccUq+nUT7geDj8MNqoTpmZCj+/BF/2Bp0Hj/4EfT6AEgW68F2YxJ6ByMnAYqCiUuotYCDwqqlRCeFAmTm5fLnxOJ3rBdCgsgXDBloba4R6ehsVER21LmjWJaPWTNnq0ONNx+zzRpoNNurWrHgFFj5i9OiTYqDNWOj2OniVMrd9USD21KKZp5SKBLphFBq7X2tdwLlTQjjf0p2niE/JZNJDFpQl0BpWvmb0sgFKVzF62y1HGmPmd+KP/0DiEXhkmXN6zu2eMpJ8+PtQvpaxYEiN9ua3KwrM3qkEZzEW/SgOlFRKhUotGlEYaK2Zse4oDSqXpkMdO6/qdJS8XFj+fxA5G1o9YQyhbJth1HxZ9x407Aetxxjj3Ld7Re3xDbB5qrHfWl3Mif967n4Vat8FVUPBy8d57YoCsWfJvv8AI4Ej2MbhkVo0opBYcyie6HOpfPhQM+eWJcjNNoZP9nwHHV8whjGUgvr3QsIR2DYTds6FfT9ApcbQajQ0fci+oY7MVFj6lLGwdPc3TH8rV1EKgjs6t01RYErrm6+Sp5Q6CDTRWmeZHUxYWJiOiIgwuxlRRGw/cYHnvt1BTq5m3d/vwqu4k67ry86ARaPg4M/QbQJ0euH6z8tKM74Atn4BZ/dAibLQYriR7CvUvvH+l78I276AUT/LEIlAKRWptQ673mP2fOL3AuUcGpEQJsrN00z+PZpB0zahNXz2cKjzkntmKnzzkJHce79/4+QORm+95UgYGw6PrYC63WHrdPgk1LiC9OAvxjBPfkfXGsM8bcdJche3ZE8PPgxj4e29QObl7Vrrfo4ORnrw4k6dTLzECwt3su34Bfo1q8p/7m9M2ZIOmrVyK+lJxkLWcRHQ/zNoPvSWL/mLlDMQ+ZUxbp9yGsoFQdjjxpWhxYobc849PGHsehkDF8DNe/D2JPh9wOfAHvLNf9dar7WzcQ8gAojTWt+0vJ0keHEnlu6M49XFe9HAf+4P4YEWgc5rPDXeuKQ//oBRdrfRHfZ/crPhwHLYOgNi1oNHCWPYJv6A0duvLtVChOFmCd6eWTSXtNaT76D98UAUINctC1OkZGTz+tJ9LN4RR2hQOT4e0oLq5Z3Yu70YC3PuN/4dNt8xFxx5eELI/cbP2f3GmPuu+cYJW0nuwk72JPhwpdT/gGVcPURzy2mSSqlAoA/wFsYKUUI4VGRMIs8v2EnchXSe716XZ+6q49wywIlH4av+kJFkVGCs0c7xbVRqBH0/NMb03XGBEmEaexJ8C9u/bfNts3ea5EfA34EbXoWhlBqDsTQgQUFBduxSCMjJzeOTPw7zyR/RVC1Xku/GtqNljTsokVsQ56KMnntuFjy6DKq2uOVL7kgxqe4tbo89V7IWaMk+pVRf4JzWOlIp1fUm+58OTAdjDL4gbYmi5WTiJcbP38H2E0k80KIaE/uHUMbbSSdSL4vbDnMfNMbGR/0CFRs4t30h7HDDBK+UelhrPVcpdd2hFa31h7fYdwegn1KqN+ANlFFKzdVaP1zwcEVRt3hHLK8t2YcCPh7SnP7Nqzk/iOMb4JvB4ONnlAko72YFzITbuFkP/vIlddcbXrllT1tr/QrwCoCtB/+iJHdRUMkZ2by2ZC9Ld54irIYfkwY3d+6J1MuiV8GC4cb0xUeWQpmqzo9BCDvdMMFrrT+33Vyltd6Q/zGlVAdToxIin23HE3l+/k7OJGfwwj31eKprbWvWU923BL4fbSxKPWIxlPJ3fgxC3AZ7TrJ+AoTase2GtNZrgDV2RyUExonUyb9HM2X1Yar5lWThk+1oWcPPmmB2zINlz0Bgaxi2wPELawhhgpuNwbcD2gMB14zDlwE8zA5MFG3pWbmM+TqC8OjzPBhajYn9Qijt7BOpl0XMgp/+BrXugiHzpPa5KDRu1oP3Anxtz8k/Dp+MseiHEKZIy8zh8a+2seVYIm8/2IQhrS2cPpt2Hn79J9TuBkO/heIlrItFiNt0szH4tcBapdSXWusYJ8YkirCUjGxGzd7GjpNJfDTYolky+W36FHIyoNc7ktxFoWPPPHhJ7sIpLl7K5pHZW9kXd5FPhragd5Mq1gaUfsGoBRPyAPjXtTYWIQrA3hWdhDDVhbQsHp65heizqUx9uCX3NKpkdUiwZTpkpUCn/7M6EiEK5JZzza43JVKmSQpHOp+aydAZm4k+l8r0R1wkuWemwObPoH5vqNzY6miEKBB7JhN/Yuc2IW7bueQMhkzfzPGENGaPbEXX+hWtDsmwbaZRQKzTi1ZHIkSByTRJYZnTF9MZNmML55Iz+GpUa9rUcvKi2DeSnQ6bpkDtuyGwpdXRCFFgMk1SWOJk4iWGfbGZpLRs5jzexroLmK5n+xxIi4fOL1kdiRB35LamSSqligG+WutkZwUo3M/x82kMm7GZ1Mwc5o5uQ7Pq5awO6U85mbDhY6jRQdY8FYWePWPw/1NKlVFKlcJYl3W/Ukq6NqJADp9LZfD0TaRn5/LtmLauldwBdn0LyXHQWcbeReFnT4JvZOux3w/8AtQERpgZlHBPB8+kMGT6ZnLzYP6YdoRULWt1SFfLzYHwD6FaS6MsgRCFnD0J3lMp5YmR4JdprbOxo1ywEPntO3WRoTM2U0zB/DFtqV/5hot8WWfvIkiKMWbOyNJ4wg3Yk+A/B45j1Idfp5SqgXGiVQi77I5NYtiMLXgXL8bCJ9tRp6Kv1SH9VV4uhH8AlRpDvXutjkYIh7hlgtdaT9ZaV9Na99aGGED+fhV2iYy5wPAZWyhTsjgLnmxHsL+LVmKMWgbnDxlXrcrap8JN2HMlayWl1Eyl1C+2+42AR02PTBR6W44m8MjMLfiXLsGCMe2sWYHJHlrDug+gQl1o1N/qaIRwGHu6Kl8CK4DLa5MdAp43KR7hJsKj4xk5exuVy3qzYExbqpYraXVIN3ZoBZzdY+u9yzV8wn3Yk+D9tdYLgTwArXUOkHurFymlvJVSW5VSu5RS+5RSE+8wVlFIfBdxklGzt1Gjgg/zx7SjYhlvq0O6Ma1h3XtQrgY0kev3hHuxp5pkmlKqAraZM0qptsBFO16XCdyttU61zcJZr5T6RWu9ueDhClemteajVdF8/Hs0ner68+nwUMpYtQqTvY6ugbgI6PsReLh4rELcJnsS/AvAMqC2UmoDEIAdpQq01hpItd31tP3I9Eo3lZWTxys/7OH77bEMahnIfx9sgqcVC2PfrnXvQ+mq0HyY1ZEI4XD2LPixXSnVBagPKOCgbS78LSmlPIBIoA7wqdZ6y3WeMwYYAxAUZOHSbKLALqZnM25uJBuPJPDCPfV49u46qMIwjzxmI8Ssh3vfltWahFu6ZYJXSnkDTwEdMXrg4UqpaVrrjFu9VmudCzRXSpUDFiulGmut917znOnAdICwsDDp4RcycUnpjJq9lWPn0/jwoWY8GBpodUj2W/c++PhDqEwKE+7JniGaOUAKf9aAHwZ8DQyytxGtdZJSajVwL0Y9G+EG9sZd5LEvt5GenctXo1rTvo6/1SHZLy4SjvwO3d8ALxedvinEHbInwTfWWjfKd3+1Umr/rV6klAoAsm3JvSRwD/BOAeMULmb1gXM8/c12/Hy8mDu6DfUquWDpgZsJ/xC8y0HY41ZHIoRp7DkLtt02cwYApVQbIMKO11XB+DLYDWwDVmqtfypYmMKVfLPlBKPnRFAroBSLn2pf+JL72X1w4CdoOw68y1gdjRCmudmKTnswxtw9gY1KqRO2+zWAA7fasdZ6N9DCQXEKF5CXp3nvt4NMXXOEu+oHMGVYKKVKFMJ128M/AC9faD3G6kiEMNXNfjv7Oi0K4fIyc3J58bvd/LjrFMPbBDGxXwjFC8M0yGudj4a9P0CH8eBT3upohDDVzVZ0inFmIMJ1JV3KYsycSLYeT+TlXg14snOtwjEN8nrWT4Li3tDuGasjEcJ0hfDva+FMJxIuMfLLrcQmpvPJ0Bbc16zqrV/kqi7EwK75xtCMb4DV0QhhOknw4oZ2nkxi9FfbyM7VzB3dhtY1C/mQxoaPjGJi7Z+1OhIhnEISvLiu3/ad4bn5OwgoXYIFo1pTO8AFF+m4HcmnYcdcaD4cylazOhohnEISvLhKbp5m6prDfLDyEE0DyzHz0TD8fd3gMv6NnxirNnV83upIhHAaSfDiitgLl3hhwS62Hk+kX7OqvDOgKSW93KA+etp5iJgFTR8Cv2CroxHCaSTBCwCW7TrFvxbvQWv48KFmPNCiWuGdKXOtTZ9CTgZ0fMHqSIRwKknwRVxKRjYTlu7jhx1xhAaV46PBLQiq4Ea1WdIvwNYZEHI/BNSzOhohnEoSfBEWcTyR5xfs5FRSOs93r8szd9VxrYuXkk/DiU3G2Hlezi1+rn2O7f75Q5CVAp1etPrdCOF0kuCLoJzcPCb/cZgpf0RTza8k341tR8saLjYFcv9SWPosZNqzeJiNKgbFiuf78TD+bTkSKjc2LVQhXJUk+CLmRMIlxi/YwY4TSTwYWo2J/UIo7UrL6mVnwG//gm1fQNVQ6P2eUfXxcrLOn7w9PP+8rzygmAv99SGEC5AEX0Rorfl+exwTlu6lWDHF5KEt6OdqV6WePwzfjYSze4xSAt0mQHEvq6MSotCSBF8EXLyUzT+X7GH57tO0rlmeSYObU61cSavDutqu+fDTC8bSecMWQr2eVkckRKEnCd7NbTqSwAsLdxKfkslLPesztkttPIq50PTHrDT4+SXYOQ9qdIAHZ8iVpkI4iCR4N5WVk8ekVYeYtvYIwRVK8f249jSrXs7qsK52Zi8sGmWU8O3yD+j8d/CQj6QQjiK/TW7oSHwq4+fvYG9cMkNaVee1vo1ca2EOrY0rS399BUqWg0eWQq0uVkclhNsx7bdeKVUdY8HuShgrQU3XWn9sVnsC0jJzmL3hGJ+uPkIJz2JMeziUextXsTqsq2VchGXPwf4lULsbPPC5lO4VwiRmdutygP/TWm9XSpUGIpVSK7XWt1ywW9yerJw85m87weTfD3M+NZN7GlXiP/0bU7mst9WhXS0uEr4bBRdjoftEaP+cTG0UwkSmJXit9WngtO12ilIqCqgGSIJ3kLw8zbJdp/hg5UFOJqbTOrg8n48Idb2LlrQ26sGsegNKV4bHfoXqra2OSgi355SBWaVUMMYC3Fuu89gYYAxAUFCQM8Ip9LTWrDkYzzu/HuDAmRQaVinD7FGN6VovwPUKhKUlwJJxEL0CGvSF/lOgpJ/VUQlRJJie4JVSvsD3wPNa6+RrH9daTwemA4SFhWmz4ynsIo4n8u6vB9l6PJGg8j58PKQ59zWtSjFXmvp4WcxGWPQ4XDoPvd6D1k+Aq30BCeHGTE3wSilPjOQ+T2v9g5ltubsDZ5J5f8VBVkWdw9+3BP/pH8LgVkF4FXfBMezUcxD+IWz9HPxqwuhVUKWZ1VEJUeSYOYtGATOBKK31h2a14+5OJl5i0qpDLN4Rh69XcV7qWZ9RHYLx8XKhaY+XpSUY655unQG5WdBiBPR8C0qUtjoyIYokM7NEB2AEsEcptdO27Z9a659NbNNtnE/NZMofh5m3JYZiSjGmUy3GdqmNXykXrM1yKRE2TYEtnxtXpjZ9yLhwqUJtqyMTokgzcxbNekAGXG9TSkY2M8KP8UX4UTJz8ngoLJDnutWlSlkXqx0DkJ4Em6fC5s8gMxlCHoSuL0NAfasjE0IgV7K6DK2NKY8Tf9xPYloWfZpU4YUe9agd4Gt1aH+VmQJbphkLWWdchIb3QddXoFKI1ZEJIfKRBO8CLqRl8erSvSzffZoWQeWYPbKV69WNAWP4Zet02DAZ0hOhXi+46xU5gSqEi5IEb7HVB8/xj0W7SUzL4qWe9Xmycy3XWjYPIDvdqB2zfhKkxUOd7nDXP6FaS6sjE0LchCR4i6Rl5vDfn6OYt+UE9Sr5MmtkKxpXK2t1WFfLzoDtXxlTHlPPQK2u0PWfENTG6siEEHaQBG+ByJgLvLBwJycSL/FEp5r8X4/6eHt6WB3Wn3KyYMfXEP4BJMcZddoHzoTgjlZHJoS4DZLgnSgrJ4+Pfz/E1DVHqFK2JN8+0Za2tSpYHdafcnNg93xY+w4knYDA1nD/Z1Czi1yBKkQhJAneSQ6eSeFvC3ay/3Qyg1oG8vp9jVxnseu8PNj3A6z5HyQchirNoc+Hxli7JHYhCi1J8CbLzdPMWn+M91YcpLR3caaPaEmPkMpWh2XQGg4sh9Vvwbn9ULERDJ4HDfpIYhfCDUiCN9HJxEv833e72HoskR6NKvHfB5vg71vC6rCMxH74d1j9JpzaARXqwICZxoVKUp9dCLchCd4EWmu+i4zl3z8ape/fG9iUgS0DXaOU77Fw+ONNOLkZygZB/0+h6RBZC1UINyS/1Q52PjWTV37Yw8r9Z2lbqzzvD2pGoJ+P1WHByW1Gj/3oGihdBfp8AC0egeIuWNtGCOEQkuAdaMW+M/zzhz2kZObwap+GPNahpvV12k/vgtX/hUO/go8/9PwvhD0Gni5Y20YI4VCS4O+Q1po/Dpxj6pojRMRcIKRqGb4d3Jx6lSwukXvuAKz5L+xfCt5lodvr0PpJKOGCtW2EEKaQBF9AObl5/LT7NFPXHOHg2RSqlSvJG/c1YlibGtYtwpGeBCe3wN7vYfdC8CoFnf8O7Z6GkuWsiUkIYRlJ8LcpIzuX7yJO8vm6o8ReSKduRV8+fKgZ9zWriqeza8gknzKWxTuxCU5shrP7AA3FS0L7Z6HD81DKhS6kEkI4lSR4O11Mz2bu5hhmbzjG+dQsWgSVY8J9IXRrUNE54+xaw/loOLERYjYZST0pxnjMsxRUb2WU7A1qC4FhRu9dCFGkSYK/hXMpGcxaf5x5m2NIycyhS70AxnWtTZua5c2d9pibDad323rntp9LCcZjPv5Qox20eRKC2kHlpjLNUQjxF2auyToL6Auc01o3Nqsds8QkpPH5uqMsiowlJzeP3k2qMLZLbXMrPl6MhR3zIGYDxG6D7EvGdr9gqNvTSOpB7YwLk1xhTr0QwqWZ2e37EpgCzDGxDYfbfyqZqWuPsHz3KYoXK8aAloE82bkWwf4mDnkkHjMWq94xD/JyoHJjaPGwkcyD2kGZKua1LYRwW2auybpOKRVs1v4dSWvNlmOJTFt7hDUH4ynl5cETnWrxeMeaVCzjbV7D8Ydg/YfGjJdiHhD6CHR8HsoFmdemEKLIsHzgVik1BhgDEBTk3MR2KimdxTviWBQZy7HzaVQo5cWLPeoxom0wZX1MrPR4dh+sex/2LYbi3tBmrDHrRXrqQggHsjzBa62nA9MBwsLCtNntpWfl8uu+03wfGceGI+fRGtrULM+4rrW5r2lVSnqZuPBG3HYjsR9cDl6+Rm+97dPgG2Bem0KIIsvyBO8MWmu2Hb/AosiT/LznDKmZOVQvX5Lx3eryYItAgiqYXCvmxGZY9x4cXmVcVdrlZWMGjE95c9sVQhRpbp3gTyZe4oftcXy/PZYTiZco5eVB7yZVGNgykFbB5c2dv641HFtnJPbj4eBTAbpNgFajwbuMee0KIYSNmdMkvwW6Av5KqVhggtZ6plntXZaWmcMve8+wKPIkm48mAtC+dgWe716XextXxsfL5O80rSF6pZHYY7eCb2WjwFfLkXLxkRDCqcycRTPUrH1fKy9Ps/lYAt9HxvHL3tNcysqlRgUf/u+eejwQWs055XqzMyB6hbFQ9eldULa6UZK3+cPgaeJMHCGEuIFCP0STlpnDvR+v42RiOqVLFKdfs6oMbBlIyxp+5i+wkXTC6K1Hr4Rja40Lk/xqQr8p0HSw1FoXQliq0Cf4UiWK06txFUKqlqFHo8rmzoLJzTZOmEb/ZiT1+Chje7ka0Hw41O0Bte+WsgFCCJfgFpnon70bmrfz5NPG7Jfo3+DIashKgWKeENwBQkcYSV1KBwghXJBbJHiHys2BuAhbL/03OLPH2F6mGjQZYCT0mp2hhMULegghxC1IggfISoMDy41l7Q7/DhlJoDyM0rvd3zCSesVG0ksXQhQqRTfBaw2ntsP2ObDne2PopVQANOgDde+BWnfJKkhCiEKt6CX4S4mwewFs/xrO7TNWPwq5H1qMMCo3FrNouT0hhHCwopHg8/Lg2BojqR/4CXKzoGoo9J0EjQcY5QOEEMLNuHeCv7yAxs65xpz1kn4Q9pjRW69c6NYgEUKI2+J+CT4nCw7+bIytH/kD0FCrq1EHpkFfuapUCFFkuE+CPxdlDMHsnm+sXVqmGnT5u3EBkl8Nq6MTQginK/wJPjMVvr7fWMO0mCc06A0tHoHadxmrJAkhRBFV+BN8CV+j/kuj+6HZECjlb3VEQgjhEgp/ggcYMMPqCIQQwuXIpG8hhHBTkuCFEMJNSYIXQgg3ZWqCV0rdq5Q6qJQ6rJR62cy2hBBCXM20BK+U8gA+BXoBjYChSqlGZrUnhBDiamb24FsDh7XWR7XWWcB8oL+J7QkhhMjHzARfDTiZ736sbdtVlFJjlFIRSqmI+Ph4E8MRQoiixfKTrFrr6VrrMK11WEBAgNXhCCGE2zDzQqc4oHq++4G2bTcUGRl5XikVU8D2/IHzBXytu5FjcTU5HleT4/EndzgWNyy2pbTWprSolCoOHAK6YST2bcAwrfU+k9qL0FqHmbHvwkaOxdXkeFxNjsef3P1YmNaD11rnKKWeAVYAHsAss5K7EEKIvzK1Fo3W+mfgZzPbEEIIcX2Wn2R1oOlWB+BC5FhcTY7H1eR4/Mmtj4VpY/BCCCGs5U49eCGEEPlIghdCCDdV6BN8USloppSqrpRarZTar5Tap5Qab9teXim1UikVbfvXz7ZdKaUm247LbqVUaL59PWp7frRS6lGr3tOdUkp5KKV2KKV+st2vqZTaYnvPC5RSXrbtJWz3D9seD863j1ds2w8qpXpa9FbumFKqnFJqkVLqgFIqSinVrqh+NpRSf7P9juxVSn2rlPIusp8NrXWh/cGYfnkEqAV4AbuARlbHZdJ7rQKE2m6XxrjGoBHwLvCybfvLwDu2272BXwAFtAW22LaXB47a/vWz3faz+v0V8Ji8AHwD/GS7vxAYYrs9DRhnu/0UMM12ewiwwHa7ke0zUwKoafsseVj9vgp4LL4CRttuewHliuJnA6McyjGgZL7PxMii+tko7D34IlPQTGt9Wmu93XY7BYjC+DD3x/jlxvbv/bbb/YE52rAZKKeUqgL0BFZqrRO11heAlcC9znsnjqGUCgT6AF/Y7ivgbmCR7SnXHovLx2gR0M32/P7AfK11ptb6GHAY4zNVqCilygKdgZkAWussrXUSRfSzgTH9u6TtYksf4DRF9LNR2BO8XQXN3I3tz8gWwBagktb6tO2hM0Al2+0bHRt3OWYfAX8H8mz3KwBJWusc2/387+vKe7Y9ftH2fHc5FjWBeGC2bcjqC6VUKYrgZ0NrHQe8D5zASOwXgUiK6GejsCf4Ikcp5Qt8DzyvtU7O/5g2/rZ0+3mvSqm+wDmtdaTVsbiI4kAoMFVr3QJIwxiSuaIIfTb8MHrfNYGqQCkK518hDlHYE/xtFzQrzJRSnhjJfZ7W+gfb5rO2P6+x/XvOtv1Gx8YdjlkHoJ9S6jjGsNzdwMcYQw2Xr87O/76uvGfb42WBBNzjWIDRu4zVWm+x3V+EkfCL4mejO3BMax2vtc4GfsD4vBTJz0ZhT/DbgLq2M+ReGCdJllkckyls44IzgSit9Yf5HloGXJ7t8CiwNN/2R2wzJtoCF21/rq8Aeiil/Gy9nR62bYWG1voVrXWg1joY4//8D631cGA1MND2tGuPxeVjNND2fG3bPsQ2k6ImUBfY6qS34TBa6zPASaVUfdumbsB+iuBnA2Nopq1Sysf2O3P5WBTJz4blZ3nv9AdjRsAhjLPc/7I6HhPfZ0eMP7F3AzttP70xxgt/B6KBVUB52/MVxpKJR4A9QFi+fT2GcdLoMDDK6vd2h8elK3/OoqmF8Ut4GPgOKGHb7m27f9j2eK18r/+X7RgdBHpZ/X7u4Dg0ByJsn48lGLNgiuRnA5gIHAD2Al9jzIQpkp8NKVUghBBuqrAP0QghhLgBSfBCCOGmJMELIYSbkgQvhBBuShK8EEK4KUnwolBSSv2slCp3G88PVkrtdXAMqY7cn22fzZVSvfPdf0Mp9aKj2xFFgyR4UShprXtro6CWu2mOcX2DEHdMErxwOUqpl5RSz9luT1JK/WG7fbdSap7t9nGllL+tZx6llJphqwH+m1KqpO05LZVSu5RSu4Cn8+3fWyk1Wym1x1ac6y7b9uVKqaa22zuUUq/bbv9bKfWEHTFvs9VXn2jbdrPYWtmeu1Mp9Z6tdrkX8G9gsG37YNvuGyml1iiljl4+LkLYQxK8cEXhQCfb7TDA11aHpxOw7jrPrwt8qrUOAZKAAbbts4FntdbNrnn+0xj1t5oAQ4GvlFLel9u1ld/Nwahhwk3aBUAp1cMWQ2uMHnhLpVRnO2J7UmvdHMjFCCgLeB2jJnlzrfUC23MbYJTybQ1MsB0LIW5JErxwRZEYSbIMkAlswkj0nTCS8LWOaa135nttsG18vpzW+nJi/jrf8zsCcwG01geAGKCebd+dMRL7cowvFh+gptb64E3i7WH72QFsx0jIdW8RW2mt9Sbb9m9usm+A5dqoS34eo2BYpVs8XwjAKDMqhEvRWmcrpY5hrMSzEaO+yl1AHYyFTq6Vme92LlCygE1vw/giOYqx2IU/8ARGYr4ZBfxPa/35VRuNuv2OiO3afcjvrbCL9OCFqwoHXsQYGgkHxgI7tJ3Fk2wnYJOUUh1tm4Zfs+/hAEqpekAQcNA2RHISGITxV0P+GG5mBfCYrVY/SqlqSqmKt4gtRSnVxrZpSL6HUzCWZBTijkmCF64qHGMd2k1a67NABtcfnrmZUcCnSqmdGL3syz4Diiml9gALgJFa68u95HCMxUTSbbcDb9Wu1vo3jGGWTbZ9LuLWSfpxYIYttlIYKwmBUda20TUnWYUoEKkmKYQFlFK+WutU2+2XgSpa6/EWhyXcjIzlCWGNPkqpVzB+B2MwzjcI4VDSgxdCCDclY/BCCOGmJMELIYSbkgQvhBBuShK8EEK4KUnwQgjhpv4fa0dj+VvXt3EAAAAASUVORK5CYII=", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXmNJREFUeJzt3QdcVeX/B/APe4MCLmQ4QFFx771npWZDTc0c9avMhi1blpZp498wy7IcDVflzD1yZ+6F4kZFQVFU9r73//o+R66gaIDAXZ/363XknnMv9z6cI/d+eKaNXq/Xg4iIiMgE2Rq7AERERER3w6BCREREJotBhYiIiEwWgwoRERGZLAYVIiIiMlkMKkRERGSyGFSIiIjIZNnDjOl0OkRHR8PDwwM2NjbGLg4REREVgEzhlpiYCD8/P9ja2lpuUJGQEhAQYOxiEBERURFERUXB39/fcoOK1KTk/KCenp7GLg4REREVQEJCgqpoyPkct9igktPcIyGFQYWIiMi8FKTbBjvTEhERkcliUCEiIiKTxaBCREREJsus+6gUVHZ2NjIzM41dDCoBDg4OsLOz47klIrJQ9pY+TvvSpUu4ceOGsYtCJahMmTKoWLEi59IhIrJAFh1UckJK+fLl4erqyg8yCwyiKSkpiI2NVfuVKlUydpGIiKiY2Vtyc09OSPHx8TF2caiEuLi4qK8SVuRasxmIiMiyWGxn2pw+KVKTQpYt5xqzHxIRkeWx2KCSg2sAWT5eYyIiy2XxQYWIiIjMF4MKERERmSwGFRPUoUMHvPzyyzC15pUlS5YYuxhERGRlGFRM0KJFi/Dhhx+q21WqVMFXX31Vaq/9wQcfoEGDBnccj4mJQc+ePUutHEREZHxnriTh4o1Uo5bBYocnmzNvb+9if86MjAw4OjoW+ftlQjUiIrJ8lxPS8NfBaCw7GI1DF+IxrHUVvP9QHaOVx9bqJgjLyCr1TV63KE0/8vXcuXN45ZVXVNNL7tEt27ZtQ9u2bdU8IgEBAXjxxReRnJxsuF9qYqRW5sknn4SnpyeeeeYZdfzNN99EjRo11JDeatWq4b333jMM6509ezbGjx+PgwcPGl5Pjt3e9NOqVSv1PLlduXJFTWe/ZcsWtZ+eno7XXnsNlStXhpubG5o3b45NmzYV+doREVHJiU/NxO+7o/DEj/+ixaQN+GhFhAopdrY2SEjNgjFZVY1KamY2ao9bU+qve3RCd7g62hepCah+/foqZDz99NOG46dPn0aPHj3w0UcfYebMmSokvPDCC2qbNWuW4XGff/45xo0bh/fff99wzMPDQ4UPPz8/HD58WD2vHHvjjTfQv39/hIeHY/Xq1Vi/fr16vJeX1x3lGjRoED799FNMnjzZEJ4WLFignlPCk5CyHD16FPPnz1fHFy9erMosrxkSElLoc0FERMUrLTMbG4/FYsmBi9h47AoysnWG+xoHlUXfBn7oVbcSfNydYExWFVTMsQlIZlqVIJG76WXSpEkqLOR0uJUP/ilTpqB9+/aYNm0anJ2d1fFOnTrh1VdfzfOc7777bp5aF6n1kDAhQUVqZ9zd3WFvb3/Ppp7HH39cvXZOrY6YO3cuBg4cqILL+fPnVWCSrxJShLyOBCA5/vHHHxfzmSIiooLI1umx43ScCidrwi8hMf1WbUmNCu7o06Ayetf3Q4C36UyWalVBxcXBTtVuGON1i5M0zRw6dAhz5swxHJPmJZ1Oh8jISNSqVUsda9KkyR3fKzUfEmqkViYpKQlZWVmqaagwypUrh27duqnXl6Air7ljxw788MMP6n6pNZElDKSJKTdpDuJyBkREpUuv1+PghXgsPXARyw/F4EpiuuE+Py9n9G5QGX0a+CG0oodJTqBpVUFFLkBRmmBMjQSM//3vf6pfyu0CAwMNt6VvSG4SJqQmRvqhdO/eXTXrSG3K//3f/xW6DPI88vrffPONqk2pW7eu2nLKJzVBe/fuvWPtHamxISKiknf6ShKWHojGsgMXcTYuxXC8jKsDHqhbSdWeNAkqC1tb0wsnuZn/p7aFk5E6UjuRW6NGjVT/j+Dg4EI91z///IOgoCC88847hmPSWfe/Xi8/ffr0UX1npDlHgop02s3RsGFD9RyyUGBO0xAREZW8S/FpWH4oWgWUwxfj89Tsd61dQdWctA0pB0d78xlLw6Bi4qQfiYykGTBgAJycnODr66tG3LRo0UJ1WB05cqSqOZHgsm7dOkydOvWuzyV9WaTfiNSiNG3aFCtWrFCdXG9/PWnKOXDgAPz9/VX/GHnd28lr9u3bV40aioiIUP1TckiTj9S4SHiR2hoJLtLhd8OGDahXrx4eeOCBYj5LRESmIz0rG1tOXDUEBBdHO7je3LTb9uq4jKgpDvEpmVh9JAZL9kfj38g45Aw0ledvF+Krak4kpLg5medHvnmW2opMmDBBNfNUr15d9fGQtkb5sN+8ebOqGZEaCzkm98uonXvp3bu3GuosAUeeSwKDBA2Z5C3HI488okYbdezYETdu3FCdX5966ql8n0/CSK9evdCuXbs8TU5Cvk9GJUln3osXL6qAJeHqwQcfLKYzQ0Rkmj5ZdRwzt0f+5+OkVkMFGEOYsTeEmrwBx17t54QdZ3XbHhnZ2Vh1+BI2Hc87Ykeac/o0rIxeYRWNPmKnONjoCzvJhwlJSEhQ/Szi4+Pv6BCalpamagaqVq1qGAVDlonXmohMRWxCGtp8uhEZWTqEVfZEZpYeKZlZSM3IVltKZrahxqM41azggT4N/fBQPdMasVOUz+/bsUaFiIiomEzfckaFlEaBZbDwuVZ3jKKRuoH0LB1SJLRkaAFGu52t5jUxHDfczkaqmjw0Wx3Lebx8zdTp0KKaz80RO4UbvWlOGFSIiIiKwdWkdPy2Uxug8GLnkHyH+soxabqRzdut6MuaWBPz6fZLRERkwn7aGom0TB3q+XuhfY1yxi6OxWBQISIiuk/XkzPwy46z6vaLnfKvTaGiYVAhIiK6TzLKR/qO1K7kic61yvN8FiMGFSIiovtceXj29pu1KZ2DWZtSzBhUiIiI7oOEFFncT4YId6t99wVdqWgYVIiIiIooMS0TM7adUbdf6BRs8uvmmCMGFbovMmutTKVfkmRa/6+++qpEX4OIqCh+2XEOCWlZqFbODb3qVuJJLAEMKlZm9uzZKFOmTLE939dff62ek4jI2iSnZ+GnrVptyuhOwcW2dg/lxQnfqEhkdWQZfidTIBMRWaM5O8/hekomgnxc1dT1VDJYo2Jipk+fDj8/P+h0txaYEn369MHw4cPV7aVLl6JRo0ZqDaNq1aph/PjxyMrKMjxWFhOUhQwrVKigHhMWFobly5dj06ZNGDZsmFpbQUKGbDkLEl6/fl2tdly2bFm4urqiZ8+eOHny5B01McuWLUPt2rXVisqyEnPupp+zZ88anjf31qFDB8PzbNu2TS2k6OLigoCAALz44otITk423B8bG4uHHnpI3S/rNM2ZM6cEzzYRUdHIFPbTt2gLD47qEAx7O36clhTrqlGRlaAyU0r/dR1cZd7kAj30sccew+jRo7Fx40Z07txZHbt27RpWr16NlStXYuvWrSpQTJkyRX3gnz59Gs8884x63Pvvv68CjoSMxMRE/Pbbb2pV5aNHj8LOzg6tWrVSfT3GjRuH48ePq+9xd3dXXyVwSDCRICILRL355ptqZWT5XgcHB/WYlJQUfPLJJ/jpp5/g4+OD8uXzzhUgwSMmJsawf+nSJXTp0kWtriykrD169FCrKs+cORNXrlxRKznLJqst55QjOjpa/fzyuhJkJLwQEZmSebvOqynzK5dxwcONKhu7OBbNuoKKhJSPjVA993Y04OhWoIdKjYYEjblz5xqCyp9//glfX1907NgR3bp1w9ixYzF06FB1n9SofPjhh3jjjTdUUFm/fj127dqFiIgI1KhRw/CYHNJUI7UcFSveGkKXE1C2b9+uwoyQmgwJHkuWLFHhSWRmZuK7775D/fr18y27hKGc55UVjaWmpWXLloZam0mTJmHQoEF4+eWX1X5ISIgKXO3bt8e0adNUDc2qVatU+Zs2baoeM2PGDNSqVavQp5yIqKTI4oE/bDmtbj/fsTocWJtSoqwrqJgJ+TB/+umnVSiQJhYJDQMGDICtrS0OHjyoAsXEiRPz9BeRYCA1HgcOHIC/v78hpBSEhBp7e3s0b97ccExqTGrWrKnuy+Ho6Ih69eoV6DmlmUpqddatW6fKLaTshw4dytOcIyuJSi1QZGQkTpw4ocrRuHFjw/2hoaHF2vmXiOh+/bEnCpcT0lHJyxmPNvbnCS1h1hVUpAlGajeM8bqFIH005AN8xYoVqmZBmnu+/PJLdV9SUpLqk9KvX787vk/6o0jfjpIiz12Q9SukaWfNmjWqZsTDw8NwXMoufWekOed2gYGBKqgQEZmyjCwdpm3SalOebV8dTvZ2xi6SxbOuoCIfsgVsgjEmCRwSRKTm4dSpU6pmQzrPCvkq/UuCg4Pz/V6p8bhw4YL60M+vVkVqRaQGJjdpWpHOuDt37jQ0/cTFxanXkY6zhbFw4UJMmDBBNeFI/5jcpOzS5+VuZZfaEynH3r17DU0/UgbpHExEZAoW7ruA6Pg0lPdwQv+mAcYujlWwrqBiZs0/Dz74II4cOYLBgwcbjktHWDkuNRCPPvqooTkoPDxc1WRIfw/pvPrII4/giy++UKHg2LFjqiZEOrLK5GlSs7FhwwbV10RG+EhfERlVJM1NP/zwg6oFkX4wlStXVscLSsogHX2lI26dOnVUZ9qccOTt7a2Ot2jRQnWeHTlyJNzc3FRwkeahqVOnqkAmZZRaF+mzIs1A0p+lJGuJiIgKKjNbh+82nVK3n2lXDc4OrE0pDRxPZaI6deqkPtylRuGJJ54wHO/evbsaarx27VpV6yAf/NIsFBQUlKdWQ+4bOHCgqhGRjrY5tShSY/Lss8+if//+KFeuHD799FN1XEbdSN8QCUHSAVaanmSUUc6In4LYs2eP6icjgalSpUqGLaeZSmp7Nm/erGp7ZMRSw4YNVfCS4dg5pByyL4FLvk9GNN0+uoiIyBiW7L+IqGup8HV3xKDmt95zqWTZ6OUTyUwlJCSoUSwyL4gMqc1NOpdKB02Zi0OaUshy8VoTUUnL1unR5YvNiLyajLE9Q1X/FCqZz+/bsUaFiIjoPyw/FK1CSllXBwxpwdqU0mQyQWXy5MmqH0XOHBtERESmQKfT45u/tb4pI9pUhZsTu3daXVDZvXu36sRZ0Dk6iIiISsuq8Es4FZsET2d7PNmqCk+8tQUVGYEiI1x+/PFHNSsrERGRadWmaOueDWtdFZ7OBR9gQBYSVEaNGoUHHnhArQnzX9LT01UHnNzbfzHjvsJUQLzGRFRS1kVcxrFLiXB3ssfw1lV5oo3AqA1t8+fPx759+1TTT0HIWjEyK2tB5F5Ij/NwWDa5xqIwQ6mJiAryR9CUDVptytBWQfBy5XuMVQWVqKgovPTSS2qyr4IOH37rrbcwZswYw77UqMjCeXdbIE/WiMlZeVcmNivI9O9kXm8iElLkGsu1lmtORFRcNh6PxZHoBLg62mFEm1uLu5KVBBWZJl0+YHKmhhcyKdmWLVvULKXSzHP7B48s0CdbQeWs5JsTVsgySUjJvRo0EVFx/CH09QZtpI8MR/Z2c+RJtbag0rlzZxw+fDjPsWHDhqn1XmSq9eL461hqUGRmVJnZNDMz876fj0yPNPewJoWIitvWk1dxMOoGnB1sMbIta1OsMqjIejJhYWF5jsnaLz4+Pnccv1/yQcYPMyIiKmzflCeaBaGcR8Fr8skCR/0QERGZkh1n4rDn3HU42tvif+1Zm2JsJjW93qZNm4xdBCIisnLf3OybMqBpACp4cq04Y2ONChER0U27z15TNSoOdjZceNBEMKgQERHdlNM35dHGAfAr48LzYgIYVIiIiADsP39djfaxs7XB8x2q85yYCAYVIiIi6Ztyc4Xkfg0rI8DblefERDCoEBGR1Tt8IR5/H4uFrQ0wqmOw1Z8PU8KgQkREVi9nheQ+DSqjiq+b1Z8PU8KgQkREVi0iJgFrj16GLAfH2hTTw6BCRERWberNvikP1K2E4PLuxi4O3YZBhYiIrNbJy4lYGR6jbr/QiX1TTBGDChERWa2pG09Brwe616mA0Iqexi4O5YNBhYiIrNKZK0n462C0uj26U4ixi0N3waBCRERW6duNp6HTA51DyyOsspexi0N3waBCRERW53xcCpYcuKhuj+7M2hRTxqBCRERW57tNp5Ct06NdjXJoEFDG2MWhe2BQISIiq3LxRioW7rugbr/UmSN9TB2DChERWZXvN51GZrYerar7oHGQt7GLQ/+BQYWIiKzGpfg0LNgdpW6/yL4pZoFBhYiIrMaUv08iI1uHZlW80aKaj7GLQwXAoEJERFZh//nrmLfrvLr9arcaxi4OFRCDChERWbysbB3eWRyuZqF9pJE/mrM2xWwwqBARkcX7ecc5HI1JgJeLA97uFWrs4lAhMKgQEZFFi4lPxRdrj6vbb/UMhY+7k7GLRIXAoEJERBZtwl9HkZyRjUaBZfB4kwBjF4cKiUGFiIgs1t/HLmNV+CXY2dpg4sN1YWtrY+wiUSExqBARkUVKzcjGuKVH1O0RbaqiViVPYxeJioBBhYiILNLUjSdx4Xoq/Lyc8RIndzNbDCpERGRxTl5OxPQtZ9TtD3rXgZuTvbGLREXEoEJERBZFr9fjnSXhaj2fLrUqoFudisYuEt0HBhUiIrIoC/ddxK7Ia3BxsMMHvWsbuzh0nxhUiIjIYlxPzsDHKyPU7Ze7hMC/rKuxi0T3iUGFiIgsxierj+FacgZqVvDA8DZVjV0cKgYMKkREZBH2nL2G+buj1O2JD4fBwY4fcZaAV5GIiMxe5s1FB8WApgFoUsXb2EWiYsKgQkREZm/mtkgcv5wIbzdHvNmDiw5aEgYVIiIyaxeup+Cr9ScNiw6WdXM0dpGoGDGoEBGRWftg2VGkZmajWVVvPNrY39jFoWLGoEJERGZr7ZFLWB9xGfay6GDfMNjYcNFBS8OgQkREZik5PQsfLNMWHXymXTWEVPAwdpGoBDCoEBGRWfp6w0lEx6fBv6wLRncKMXZxqIQwqBARkdmJiEnAjG2R6vaHfcLg4mhn7CJRCWFQISIis6LT6fHO4sPI1unRM6wiOoaWN3aRqAQxqBARkVlZsCcK+87fgJujHcY9xEUHLR2DChERmY2rSemYvOqYuj2mW01U8nIxdpGohDGoEBGR2ZCVkeNTM1G7kieGtgwydnGoFDCoEBGRWdhxOg6L9l2ETJUiiw7ac9FBq8CgQkREJi8jS4d3lxxWt59oFoiGgWWNXSQqJQwqRERk8n7ceganryTD190Rb3TnooPWhEGFiIhM2vm4FEzZoC06+O4DteHl6mDsIlEpYlAhIiKTpdfr8d7ScKRn6dCqug/6NPAzdpGolDGoEBGRyVoVfgmbT1yBo50tPuSig1aJQYWIiExSYlomxv+lLTr4bIfqqF7O3dhFIiNgUCEiIpP0xboTuJyQjiAfVzzfobqxi0NGwqBCREQmJ/xiPH7+56xh0UFnBy46aK0YVIiIyKRk31x0UKcHHqrvh3Y1yhm7SGREDCpERGRS5u48h4MX4uHhZI/3Hqhl7OKQkTGoEBGRyYhNTMOnq4+r26/3qInyns7GLhIZGYMKERGZjI+WRyAxPQv1/L0wqDkXHSQGFSIiMhGrDsdg2cFo2Mqig33rwk5ukNWzt/ozQERERnUpPg0TV0bgr4PRan9IiyDU9ffiVSGFQYWIiIwiM1uHWdsj8fX6k0jOyIaNDTCoeSDe6sUOtHQLgwoREZW6f05dxbhlR3AqNkntNwwso+ZLCavMmhTKi0GFiIhKTUx8Kj5aEYEVh2LUvrebI8b2DMWjjfxhyz4plA8GFSIiKnEZWTrM3B6JKRtOIiUjW3WYHdwiCK92rQkvVwdeAborBhUiIipR205exfvLwnH6SrLabxxUFhP61EEdPzbz0H9jUCEiohIRfSMVE6WZ57DWzONzs5nnETbzUCEwqBARUbE388zYpjXzpGZqzTxPtqyCV7rWgJcLm3mocBhUiIio2Gw9eQXvLzuCMzebeZqoZp4w1Pbz5FmmImFQISKi+3bxRio+Wn4Uq8IvqX1fdye81TMU/RpVho1MkEJkjmv9TJs2DfXq1YOnp6faWrZsiVWrVhmzSEREVAjpWdn4duMpdPm/zSqkyLT3w1pXwd+vtccjjf0ZUsi8a1T8/f0xefJkhISEQK/X4+eff0afPn2wf/9+1KlTx5hFIyKi/7D5xBV8sOwIIq9qzTxNq2jNPLUqsZmHio+NXhKCCfH29sZnn32GESNG3HFfenq62nIkJCQgICAA8fHxqkaGiIhKp5nnw7+OYvURrZmnnIcT3u4Vir4N2MxDBSOf315eXgX6/DaZPirZ2dn4448/kJycrJqA8jNp0iSMHz++1MtGRERaM8+PW85g6sZTSMvUqWaep1pVwctdQuDhzNE8ZKE1KocPH1bBJC0tDe7u7pg7dy569eqV72NZo0JEZBw7TsfhrUWHcDYuRe03q+qt1uapWdGDl4Qsu0alZs2aOHDggCrsn3/+iaFDh2Lz5s2oXbv2HY91cnJSGxERlQ6dTq9qUL5afwI6PVDewwnvPFALvev7saMsWUeNyu26dOmC6tWr44cffijWREZERIUTl5SOlxccwNaTV9X+Y439Me6h2mzmIeuqUbmdTqfL02GWiIhK3+6z1zB67n5cSkiDs4MtPupbF4829ueloFJn1KDy1ltvoWfPnggMDERiYqLqn7Jp0yasWbPGmMUiIrLqpp4ft57Bp2uOI1unR/VybvhuUGP2RSHrDCqxsbF48sknERMTo6qAZPI3CSldu3Y1ZrGIiKzSjZQMvPbHQayPiFX7fRr44eOH68LNyeQq38mKGPV/34wZM4z58kREdNOBqBsYNWefmiPF0d4WHzxUBwObBbDDLBkdYzIRkRWT8RSz/zmLj1dGIDNbjyAfV3z7RCOEVfYydtGIFAYVIiIrlZCWiTf/PGRYSLBnWEV88mg9eHLyNjIhDCpERFYo/GI8Rs3dh3NxKXCws8E7vWphaKsqbOoh8w8qUVFR6j+yLCgodu3apUbryARtzzzzTEmUkYiIirGpZ96uKHzw1xFkZOlQuYwLvh3UCA0CyvAck0myLew3PPHEE9i4caO6fenSJTVCR8LKO++8gwkTJpREGYmIqBgkp2fhlQUH8PbiwyqkdKlVHitebMOQQpYVVMLDw9GsWTN1+/fff0dYWBj++ecfzJkzB7Nnzy6JMhIR0X06cTkRvaduw5ID0WoxQVnt+Mcnm6CMqyPPLVlW009mZqZhvZ3169ejd+/e6nZoaKiaD4WIiEzLn3sv4N0lh9WKxxU9nTH1iYZoUsXb2MUiKpkalTp16uD777/H1q1bsW7dOvTo0UMdj46Oho+PT2GfjoiISkhqRjbe+POgmsRNQkrbEF/V1MOQQhZdo/LJJ5/g4YcfxmeffaZWOq5fv746vmzZMkOTEBERGdfpK0lqArdjlxJhawO80qUGRnUMhq3sEFn66snZ2dlq5cOyZcsajp09exaurq4oX748SgtXTyYiutNfB6MxduEhJGdkw9fdCVMGNECrYF+eKrKe1ZPt7OzyhBRRpUqVojwVEREVk/SsbHy0PAK//ntO7beo5o0pAxuivIczzzGZLU74RkRkAc7HpagJ3A5fjFf7ozsF46XOIbC3K3RXRCKTwqBCRGTGktKzsHj/RXy6+hgS07JQ1tUBX/ZvgA41S68ZnqgkMagQEZnpFPhzdp7HsgMXVV8U0SSoLL55oiEqebkYu3hExYZBhYjITKRkZKmOsnN3nsfBC1oTj6hWzg2DmwdhSMsgOLCphyxMkYLKhg0b1BYbGwudTpfnvpkzZxZX2YiICEBETIIKJ0v2X0RiepY6J7KQYI+wSniiWaDqNCtrsBFZokIHlfHjx6s1fZo0aYJKlSrxl4OIqIQma1t+KBpzd53H/vM3DMer+LhiYLNAPNrYHz7u2izhRJas0EFFZqWVNX2GDBlSMiUiIrJiJy8nqr4ni/ZdQEKaVntib2uD7nUq4onmgWhZzYeTtpFVKXRQycjIQKtWrUqmNEREVigtMxurwmNU887us9cNx/3Luqjak8ea+HMuFLJahQ4qI0eOxNy5c/Hee++VTImIiKzEqdgkzNt1Hgv3XcCNlEx1TFY27lKrPJ5oHoS2wb6sPSGrV+igkpaWhunTp6uVk+vVqwcHB4c893/xxRdWf1KJiO41e+yaI5cx599z2Bl5zXDcz8sZA5oFon/TAFTw5EyyREUOKocOHUKDBg3U7fDw8Dz3sdc5EVH+zl5NVrUnf+y9gGvJGeqYrA/YKVRqTwLRvkZ5VZtCRPcZVDZu3FjYbyEislp7z13DF+tOYPupOMOxCp5O6N80EAOaBsCvDCdnI7oXTvhGRFRCDkbdwKCfdiItUweZ5qR9jXJq3hOpReEaPETFGFT69eunhiTLUsxy+14WLVpUwJcmIrJcUddSMOLn3SqktAn2xaR+dRHg7WrsYhFZZlDx8vIy9D+R20REdHc3UjIwdNYuXE3KQO1Knvh+SGO4O7ECm6gobPR6vR5mKiEhQQWn+Ph4VdtDRGQKo3qGzNiFXZHXUMnLGUtGteYoHqL7+Py2vee9RERUYDqdHq/9cUiFFA8ne8wa1pQhheg+MagQERWTz9YeV6sby5T30twTWpE1vUT3i0GFiKgYzNl5DtM2nVa3Jz9SD62DfXleiYoBgwoR0X3aeCwW7y3RJsB8uUuIWtmYiEwgqMh0+kRE1iz8YjxGzd0HnR4qoLzUOcTYRSKy7qCi0+nw4YcfonLlynB3d8eZM2fUcVmkcMaMGSVRRiIik3ThegqGzd6NlIxsNVfKxw/X5VIiRMYOKh999JGa/O3TTz+Fo6Oj4XhYWBh++umn4i4fEZFJik/NxLBZu3ElMR2hFT3w3eBGcLRnazpRcSv0b9Uvv/yiVk8eNGgQ7OzsDMfr16+PY8eOFXf5iIhMTkaWDs/+uhcnY5PUuj0zn2oKT+e8K8kTkZGCysWLFxEcHJxvk1BmZmYxFYuIyDTJHJljFx7CjjNxcHO0UyGFCwsSmVBQqV27NrZu3XrH8T///BMNGzYsrnIREZmkL9edwKL9F2Fna4PvBjdGHT8uK0JUkgq9+MS4ceMwdOhQVbMitSiyCOHx48dVk9Dy5ctLppRERCbg991RmPL3KXV7Yt8wtRoyEZlYjUqfPn3w119/Yf369XBzc1PBJSIiQh3r2rVryZSSiMjItpy4grcWH1a3X+gYjAHNAo1dJCKrUKTlPNu2bYt169YVf2mIiEzQ0egEPD9nH7J1evRt4IdXu9UwdpGIrMZ9rTuelJSkmn9y4yrGRGRJYuJTMXz2biSlZ6FFNW988mg9zpVCZMpNP5GRkXjggQdUs48s0Vy2bFm1lSlTRn0lIrIUiWnaXCmXEtIQXN4dPwxuAif7W9MyEJEJ1qgMHjxYDc+bOXMmKlSowL8siMgiZWbrVHPPsUuJ8HV3wqynmsLLlXOlEJl8UDl48CD27t2LmjVrlkyJiIiMTP4Ye3vRYWw9eRUuDnYqpAR4uxq7WERWqdBNP02bNkVUVFTJlIaIyAR88/cp/LH3AmxtgKlPNERdf86VQmQ2NSqyns+zzz6r5lGR9X0cHPJWhdarV684y0dEVurwhXicvpKE+gFlUMXHtdSamRfuvYAv1p1Qtyf0CUPnWhVK5XWJqJiCypUrV3D69GkMGzbMcEzeQKSqVL5mZ2cX9imJiPJYHR6D0fP2IzNbr/a93RzRKLAMGgaWRaPAsqgf4AVXx/satJivf05dxZsLD6nb/2tfDYNbBPHKEBlZoX/Thw8frqbKnzdvHjvTElGxW7z/Al7745CasyTIxxUx8Wm4lpyB9RGxahMyfX2tSh4qtDQO0sKLf1mX+6p1OX4pEf/7bS+ydHo8WK8S3uweWow/FREVlY1eqkIKQYYlS4fa/BYmLG0JCQlqiHR8fDznbyGyAHN2nsO7S8Ih70qPNvbHJ4/UQ5ZOpyZc23vuOvafv4F956+r8HI7GZkjtS4quASVRd3KXnB2KNhQ4ssJaXj42+2Ijk9D0ypl8euI5gX+XiKLlXodWDgS6DYRKB9qtM/vQteodOrUyWSCChFZjp+2nsFHKyLU7aEtg/D+Q3Vga2sDO1s71eQjW47oG6kqsOw7dwN7z1/H0eh4XE1Kx9qjl9UmHOxsUNvPS4WXnJqX/FY5Tk7PUhO6SUip5uuG6UOaMKQQZWcCvw8FIjcDCdHAs9sB20KPvykWhQ4qDz30EF555RUcPnwYdevWvaMzbe/evYuzfERk4aRSd8qGU/hyvdaB9bkO1fFG95r3bMaRwCHbg/X81H5aZjbCL8arWhcVYM7fwJXEdByMuqG2WdvPqsdV9HRGoyAtuEitS2hFD4yauw9HohPg4+aI2cOaoaybYyn95EQmSq8HVr6mhRRHd6Dfj0YLKUVq+rG9R2FLuzMtm36IzJu8/UxedQw/bDmj9l/rVgMvdAoplue9cD2n1kULLkdjElS/l9xk+LEccnawxfxnWqJBQJn7fm0is7fjW2DN2/KpDgycD9TsUewvUaJNP7ev7UNEVBQ6nR7jloXjt3/Pq/33HqyNEW2qFsvJlD+aZII22fo0qKyOpWRk4dCF+DzhRTrpSliZMqAhQwqROL4KWPOOuonuE0skpBRW8Y/vIyL6D1nZOryx8BAW7bsIaeGZ9HBdDGgWWKLnTYYzt6jmo7acWpdzcSmQOpaqvm68ZkSXDgN/jpDfDqDxU0CL503inBQoqEyZMgXPPPMMnJ2d1e17efHFF4urbERkgTKydHh5wX6sPHxJDTP+4vH6hlqP0iS1LlUYUIg0iZeBuQOAzGSganug1+fySwJTUKA+KlWrVsWePXvg4+Ojbt/1yWxscOaM1tZcGthHhci8SKfX537bi43Hr8DRzhbfPNEQ3etUNHaxiKxbZiow+wHg4l7AJwQYuQ5wuTXKriQUex+VyMjIfG8TERVUUnoWnv55D3aciVOdV38Y0gTta5TjCSQyJp0OWPysFlIknDyxoMRDSmEVerzRhAkTkJKScsfx1NRUdR8R0e3iUzMxZMZOFVLcnezx87BmDClEpmDTJODoEsDWAej/G+BTHaam0MOT7ezsEBMTg/Lly+c5HhcXp45xeDIR5XlvSErHkBm71PBgLxcH/DK8mVpokIiM7NDvwKKntdt9vgMaDiq1ly7R4ck5iw/eTmar9fb2LuzTEZEFuxSfhsEzduJUbBJ83R3V1PS1Kt37TYmISsH5f4Glo7TbbV4p1ZBSWAUOKmXLllUBRbYaNWrkCStSi5KUlIRnn322pMpJRGYm6loKBv20E+evpaCSlzN+G9kc1cu5G7tYRHT9LDB/EJCdAYQ+CHQaZ9LnpMBB5auvvlK1KbJ68vjx41WVTQ5HR0dUqVIFLVu2LKlyEpEZOX0lCYN/2qkWDwz0dsWckc3V5GtEFrdon3zoy3Yt8ubtm1/TEoAmw4AObwH2TjAZafHA3P5AylWgUn2g33SjTo9fIn1UNm/ejNatW8Pe3vhzxXF4MpHpiYhJUB1nryZlILi8uwopFTydjV0sosLTZQPxF26FkZwQkhNK0m7893OUrwP0+wGoWNf4VyA7C5j7OHB6A+BRCXj6b8BTWy+rtJVoH5X27dvfT9mIyIIdiLqBoTN3qVE+tSt54tcRzeDjbkJ/TRLdLj0RuH7uzhAi+zeiAF3mvc+ZW3mgbBXAu6r2tezNr4nRwMo3gNgjwPSOQIexQOuXATsj/pG/5i0tpDi4amv4GCmkFJbxq0WIyCLsPBOHET/vUfOlNAwso1YillE+RCYl9hiw/Wsg7qQWSqQJ5F5k2G7ZoFsBJE8oqQI43mP5hSrtgOUvA8eWA39/qK2j8/APgG8wSt3O6cCu6dpCg9Lc49cA5qLQTT+mhE0/RKZh84kr+N+ve5CWqUPLaj74aWgTuDnx7yAyMbER2gysKXF5j7t4518rIpvUOtjaFf015SP24Hxg1RtAegJg7wJ0HQ80fbr0+oacXA/MfQzQ64Au44E2L8OcPr8ZVIjovqw5cgmj5+5HRrYOHWuWw7TBjeHscB9v7EQlIe40MKsnkHQZ8GuoDcnNCSPOtwaHlBjp6yLDgc9s0varttPmLikTULKve/koMKMbkJEINBgM9JlqEmv4FCaoFDrOyaifxMTEO44nJyer+wpj0qRJaNq0KTw8PNRkcX379sXx48cLWyQiMpKlBy7i+Tn7VEjpVbeimhafIYVMjvRB+bm3FlIqhAGDFwG1+2ijXkojpAgvf2DwYm2xP6lVidwCTGsF7J+j1bqUhKQrwLz+WkgJag08+KVJhJTCKnRQ+fnnn9V0+beTY7/88gsKO4Jo1KhR+Pfff7Fu3TpkZmaiW7duKvQQkWmbt+s8Xl5wANk6Pfo1qowpAxrC0d60hzmSFUqIBn7pDSRcAHxrAEOWAK5GmpxUmnqaPQ08tx3wb6Y1BS19Hpj/BJAUW7yvlZkGLBgE3DivNWXJ9Pj2jjBH9oWpppHuLLJJjYqzs3OeCd9Wrlx5x7T6/2X16tV59mfPnq2eY+/evWjXrl2hnouISoe8B0zfcgaTVh1T+4NbBGJC7zDY2prfX2pk4aRG4Zc+2ige+bB+chngbgILYcp6OsNXa516N34MHF8JRO3Uajykpud+6fXAstHac0qN0aA/jBfOSjOolClTJs/MtLeT4zIR3P2Qtipxt6n409PT1ZY7PBFR6Ym+kYo3Fx7C1pPaSIln2lXDWz1D811Wg8ioUq4Bv/YFrp4APP2BocsAz0qmc1Gkg27bMUBIN2314suHgd+fBOo+DvT69P5WMN7yGXD4d8DWHnj8F8A3BOaswJ1ppZlGHtqpUycsXLgwT5iQmWmDgoLg51f0Mdk6nQ69e/fGjRs3sG3btnwf88EHH+QbhgrSGYeIik5+9//YewEf/nUUielZcLK3VQFlaKsqDClkemT2ValJid4PuFcAhq0yyVWBDbIygM2fANu+0EbmyGRs0uk1uEvhnyt8IfDnzf6iD30NNH4KpqhER/2cO3cOgYGBxf7m9Nxzz2HVqlUqpPj7+xe4RiUgIIBBhagExSak4a1Fh7HhmNaGLnOkfP5Yfa7bQ6YpIxn4tR8Q9S/g6gM8tRIoHwqzcGEPsPh/QNwpbb/JcKDrh4CTe8G/X4ZfZ6UBLV8Auk+EqSrRUT8RERHYvn27Yf/bb79FgwYN8MQTT+D69etFKvALL7yA5cuXY+PGjXcNKcLJyUn9QLk3IioZ8jeMjOrp+uUWFVIc7WzxZo9Q/PlsK4YUMk3SgXTeQC2kSN+MIYvNJ6QI/ybA/7YCzW8u8LtnJvB9a+Dcjv/+XplFV352CSk1egBdJ8BSFDqovP7664a+IYcPH8aYMWPQq1cvREZGqtuFfSOUkLJ48WL8/fffqFq1amGLQ0Ql4GpSOp77bR9emn9ATYcfVtkTf41ug+c6VIcdO82SqTafSB+PyM2Ao7s2BFmGH5sbR1eg5ydax1+vAK0jsMz/svY9LYjdbRmAeQOA5Fht+PUjP93fJHUmptBTR0ogqV27trotfVUeeughfPzxx9i3b58KLIUhQ5Pnzp2LpUuXqrlULl26pI5LdZCLi0thi0ZExWDV4Ri8uyQccckZsLe1wehOIXi+Y3U42HHoMZkoWWxv4Qjg5BptjpInFmi1E+asWnttGPPqt4EDvwH/TAFOrgMe/j7v9PeycOKfI4DL4Vp/HPnZnTxgSQr9ziMdZ1NSUtTt9evXq3lPhHSuLewonGnTpqn2qQ4dOqBSpUqGbcGCBYUtFhHdpxspGXhx3n48N2efCimhFT2wZFRrvNQlhCGFTJdOp81FErEMsHMEBswBqrSBRZDmq77fAgPmaYsfXokAfuoMbP5UC2dCalpUQHPWHicTy1mYQteotGnTRjXxtG7dGrt27TKEihMnTtyzf0l+zHiZISKLsiHiMsYuOowriemQlh1p4nmxcwic7C2n+pgskHyGyKJ/hxbcGoob3BkWJ7QXENAcWPEKcHQpsHGitsBhSFfg32+1x/SdBvg3hiUqdI3K1KlTYW9vjz///FPViFSuXFkdlxE7PXr0KIkyElEJSUjLxGt/HFSrHktIqV7ODYueb43Xu4cypJDph5TVY4F9PwM2tkC/H4GaPWGx3HyAx34G+v2k1bRE79OGNIuO7wJh/WCpuCghkZXacuKKmrwtJj5NLf8xsk1VvNqtJtfqIfMIKRvGA9u+vFWb0OAJWNWyAEtfAE5vAOoP1H5+M5t0sTDDk4u0Dvvp06cxa9Ys9fXrr79W095LjYrMr1KnTp2ilpuISkFSehY+XhmBuTvPq/0gH1c1L0rTKuY7xTZZmS2f3wopD3xhXSFFePoBgxdq6/iUCTS7kFLiTT8yQ23dunWxc+dOLFq0CElJSer4wYMH8f7775dEGYmomPx7Jg49v95iCClDWwZh1UttGVLIfPwzFdj4kXa7+8dA0xGwSjY2QNkgiw8pRQoqY8eOxUcffaRWO5YRQDlkan1ZBZmITE9qRjbG/3UEA6b/i6hrqahcxgVzRzbH+D5hcHUsUsUqUenb/ROw9p1b/TJajuJVsAKFfoeSSd5k7pPbSfPP1avaQmVEZDr2nruuOsxGXk1W+wObBeDtXrXg4exg7KIRFdyBucCKV7XbbcYA7V7j2bMShQ4qsopyTEzMHbPI7t+/3zACiIiMLy0zG1+uP4Eft5yBTg9U8HTC5EfqoWPN8sYuGlHhF9pberP2pPlzQOdxVtHkQUUMKgMGDMCbb76JP/74Qy1MKKsey9o/r732Gp588snCPh0RlYDDF+Ix5vcDOBmr9SHr16gy3n+wDrxcWYtCZubYSmDRM9qqwrIScI9JDClWptBBRabLl6nvZdXi7OxsNZ2+fJVFCd99992SKSURFXhEz7RNp/D95jPI1unh6+6Ejx8OQ7c6FXkGyfyc2gD8MRTQZQH1+gMPfMmQYoWKPI/K+fPnER4erkb9NGzYECEhITDlcdhEliwjS4d5u85jyoaTavp78WC9SpjQJwzebrc6vROZjbPbgN8eBbJSgVq9gUdnAXbs+G0pSnweFSFzpkitipAmICIqfTqdHssPx+DzNcdx/pq2BldVXze82SMUPcJYi0JmKmo3MLe/FlJCugOPzGBIsWJFWg51xowZCAsLg7Ozs9rk9k8//VT8pSOiu9p28ip6f7tNLSQoIUWaeT7qG4a1r7RjSCHzFXMQmPMIkJEEVG2vrd9jz1pBa1boGpVx48bhiy++wOjRo9GyZUt1bMeOHXjllVdUc9CECRNKopxEdFP4xXh8svoYtp7UpgNwd7LH/9pVw/A2VeHmxKpxMmOXjwC/PgykxQOBLYGB8wAHZ2OXisytj0q5cuUwZcoUDBw4MM/xefPmqfBSmnOpsI8KWZPzcSn4fO1xLDsYrfYd7GwwuEUQXugYDB93J2MXj6h4RvdkJAJ+jYAnlwLO7HtoqUq0j0pmZiaaNGlyx/HGjRsjKyursE9HRP/halI6pv59CnN2nkNmtvZ3Rd8GfmoBwQBvV54/Mm86HbD1c2DjRG0/qA3Q/1eGFCp6UBkyZAimTZummn9ymz59OgYNGlTYpyOiu0hOz8JPWyMxfctpJGdkq2PtapTDG91rIqyyF88bmb/0JGDJc0DEMm2/2TPa+j12nO+HChlUxowZY7gtI3yk4+zatWvRokULdUwWKJT+KZzwjej+ZWbrMH/XeXy94ZSqTRF1K3thbM9QtA725Skmy3AtEpg/CIg9Atg6AA9+ATTipKFUxKAi0+Pf3swjTp8+rb76+vqq7ciRIwV5OiLKh3QXW3FzqPHZOG2ocZCPK17vXhO9wirB1pbTAJCFOLMJ+OMpIPU64F4B6P8bENDM2KUicw4qGzduLPmSEFmxf05dxeTVx3DoQrza93V3xEudQzCgWSAc7Io0iwCR6ZGxG/9OA9a+C+iztU6zA+YAnn7GLhmZMI5lJDKiI9Ey1Pg4tpy4ovbdHO3wTLvqGNmWQ43JwmSmActfAQ7O1fbrPwE8+CWHH9N/YlAhMoKoayn4v7XHseTAraHGg5oH4YVOwWriNiKLkhANLBgMXNwL2NgB3T4CWjzHdXuoQBhUiEpRnAw13ngKv/17a6hx7/oy1LgGgnzceC3I8kTt0kJK0mXApSzw2GygWgdjl4rMCIMKUSnZe+4aRv68B9dTMtV+2xBftSYPhxqTxdr3C7DiVSA7AyhfGxgwF/CuauxSkZlhUCEqBWuPXMLoefuRnqVDaEUPvPtAbbQJ4VBjslDZmcCat4Fd07V9Wf247zTAyd3YJSMzxKBCVMLm7jyPd5cchk4PdA4tj2+eaAhXR/7qkYVKvqoNPT67Vdvv+A7Q9jXAlqPXqGj4bklUgvOifLn+JKZsOKn2+zcJwMSHw2DP4cZkqWIOaZO4xZ8HHD2AftOB0F7GLhWZOQYVohKQla3Du0vCMX93lNp/sVMwXulaQ83sTGSRwhcCS0YBWamAdzVgwDygfKixS0UWgEGFqJilZmRj9Lx9WB8RC5lM9sO+YWroMZFF0mUDf38EbLu5/lv1zsCjM7QRPkTFgEGFqBhdS87AiJ93Y//5G3Cyt8WUgQ3RvU5FnmOyTGnxwMKRwMm12n6rF4EuHwC2dsYuGVkQBhWiYpzEbeisXThzJRleLg6YMbQJmlTx5vkly3TlBDB/IBB3CrB3BnpPBeo9ZuxSkQViUCEqBkejE1RIuZKYDj8vZ/w8vBlCKnjw3JJlOrFGq0lJTwA8/YEBvwF+DY1dKrJQDCpExbCg4P9+3YvE9CzUrOChQkpFL2eeV7LMRQW3/p/WJwV6ILAV8PgvgHs5Y5eMLBiDCtF9WHYwGq/+fkBNh9+sqjd+fLKJavYhsjgZycDSUcCRxdp+kxFAj8mAvaOxS0YWjkGFqIhmbIvEh8uPqtu96lbEF483gLMDOxGSBTq7DVg2Grh2BrB1AHp9BjQZZuxSkZVgUCEqJJ1Oj8mrj2H6ljNqf2jLIIx7qA7sZCwykaWN6ln3PrB3lrbvUQl4dBYQ1NLYJSMrwqBCVAgZWTq88edBLDkQrfbf6FETz7WvzoncyPIcXwUsHwMkav/X0XgY0HU84Oxl7JKRlWFQISqgpPQsPPvrXmw7dRX2tjb45JF6eKSxP88fWZakK8DqN7WZZoXMMvvQFKBqW2OXjKwUgwpRAcQmpmH47N0Iv5gAV0c7fDeoETrULM9zR5Y1oufQ78DqsUDqNcDGFmg1GujwFuDgYuzSkRVjUCH6D5FXk/HkzJ2IupYKHzdHzHyqKeoHlOF5I8txIwpY/gpwap22X6Eu0Ocbzo1CJoFBhegeDkTdUDUpMjV+oLcrfhneDFV83XjOyDLodMCeGcD6D4CMJMDOEWj/JtD6JcCOw+zJNDCoEN3FxmOxeH7OPqRmZqNuZS9Vk1LOw4nniyxnCvy/XgTO79D2A1oAvb8BytUwdsmI8mBQIcrHH3uiMHbRYWTr9Ggb4otpgxvD3Ym/LmQBsjOB7V8Dmz8BsjMAR3dtIUGZwM3W1tilI7oD33mJctHr9fh24yl8vvaE2u/XsDImP1IPjvZ8AycLEL0fWDoauHxY2w/uAjz4FVAmwNglI7orBhWim6T25INlR/Drv+fU/rPtq+PNHjU5RwqZv8xUYNMk4J+pgD4bcCkL9PgEqPc4YMOJCsm0MagQAbienIE3Fx7C2qOX1fv2uAdrY1jrqjw3ZFnT34uwR7SQwoUEyUwwqBCsvaln4b6L+HhlhBrZ42hniy/618eD9fyMXTSi4p/+/oEvgNBePLNkVhhUyGqdik3Cu0sO498z19R+zQoemPRIXTQKLGvsohHdH05/TxaEQYWsTlpmNr7bdBrfbzqNjGwdnB1s8VLnGhjZtioc7Nhplixs+nsZclyljbFLRlRkDCpkVbafuop3l4Sr2WZFh5rl8GGfMAR4uxq7aERFx+nvyYIxqJBVuJqUjokrIrB4/0W1X97DCe8/VAe96lbkqB4yb9fPAiteBU6t1/Y5/T1ZGAYVsmg6nR6/74nCpFXHEJ+aqUb0PNkiCK92rwlPZ04RTmY+cduOqcCmT4CsVE5/TxaLQYUs1onLiXh70WHsOXdd7deu5ImP+9VFAy4oSOYuahfw18tA7BFtv0pbbeI232Bjl4yo2DGokMVJzcjGlL9P4sctZ5Cl08PV0Q5jutbAU62qwJ6dZcmcpd4ANowH9siQYz3g4g10/xioP4ATt5HFYlAhi7LpeCzeWxqOqGupar9r7Qr4oHcdVC7jYuyiEd1fZ9kji4DVbwFJl7VjDQYDXScAbj48s2TRGFTIIsQmpGH88qNYcShG7VfyclYBpXudisYuGlHxdpb1CQEe/BKo2pZnlqwCgwqZ/fo8c3eew6erjyMxPQu2NlBT37/StQZXOyYL6Cz7LbBp8q3Osm1fA9q8DNg7Gbt0RKWGQYXM1pHoeLy9OBwHo26o/Xr+Xvj44boIq+xl7KIRlUBn2S8B3xCeWbI6DCpkdpLTs/DV+hOYuf2sqlFxd7LH691rYnCLINhJlQqRWXeWnQDsmZmrs+xEoP5AdpYlq8WgQmZl/dHLGLc0HNHxaWpfJmwb92AdVPRyNnbRiO6zs+xiYPXYXJ1lBwFdP2RnWbJ6DCpkFmLiU/HBsiNYc0R7E5dRPB/1DUPH0PLGLhpRMXSWfQ04tU7b9wnW5kRhZ1kihUGFTJper8e8XVH4eGUEktKzVNOOLB74UucQuDryvy9ZWmfZV4E2r7CzLFEufKcnk3XxRirGLjyErSevqv2GgWVUZ9lalTyNXTSi+xO1G1j+MnA5XNtnZ1miu2JQIZOsRZm/O0otIii1KE72tqqzrAw7ZmdZMmtp8cB6mVmWnWWJCopBhUxKtNSiLDqMLSeuqP1GgWXw2WP1Ub2cu7GLRnR/nWWPLgFWvcnOskSFxKBCJlOLIqscf7Q8Qk3cJrUor3WrieFtWItCZu7aGS2gnFybq7OszCzbztglIzILDCpkEiN6xi48jM03a1GkL8rnrEUhc3flBLDtC+DQ74A+W+ss22aM1lnWgcPpiQqKQYWMWovyx54L+HD5UVWL4qhqUWpgRJtq7ItC5ivmELD1/4CjS7VJ20RwF6D7JKBcDWOXjsjsMKiQUVyKT8PYRYew6bhWi9IgQKtFCS7PvihkxiN5tn4OnFh961jog0DbMUDlxsYsGZFZM2pQ2bJlCz777DPs3bsXMTExWLx4Mfr27WvMIlEp1KL8ufcCJkgtSppWizKmaw083Za1KJRLVgaweTIQdxqo9RBQowfg5G6anWTPbgO2fAZEbtaO2dgCdR7W5kSpUMfYJSQye0YNKsnJyahfvz6GDx+Ofv36GbMoVEq1KG8vPoy/j8Wq/fpSi/JoPYRU8OD5p1sSLwN/DAXO79D2ZbSMvQtQoxtQpx8Q0g1wdDV+QDm1XgsoUTu1Y7b2QL0BWh8U32Djlo/Ighg1qPTs2VNtBZWenq62HAkJCSVUMiruWpSF+y5iwl9HkCC1KHa2eLlrCJ5pWw32drY82ZR31eAFQ4CkS4CTp7YYn0wtLyNnpM+HbA5uQM2eQFg/re+HvVPpnUGdDji2XGviiTmoHbNzAhoNAVq/BJQJ5NUksuY+KpMmTcL48eONXQwqhMsJaXh70WFsuFmLUs/fS/VFqcFaFLq9hkImQZNhvLpMoFwoMGAu4FNdu09CwZFFQPhiIP48EP6ntkmYCX1Aq2mp1gGwdyyZ85qdpb2+dJK9ckw75uAKNBkOtBoNeFTk9SQqITZ6+XPXBNjY2PxnH5X8alQCAgIQHx8PT09Oq25K5L/V4v0X1UKCObUoL3UJwf/asRaFbpOZBqx4FTjwm7Zfuw/Q57v8+6TI29XFvUD4Im214cToW/c5l9H6s0hNS5V2gJ198fSVOTQf2PoFcD1SOybhqNkzQIvnubIxURHJ57eXl1eBPr/NqkbFyclJbWTaYqUWZXE41kdoKx3XrazVotSsyL4odJsbUcDvQ4Do/Von1C4fAK1elL9c8j9Vcty/ibZ1+0jrHxK+UGsSSo4F9v+qba6+QO3eWk1LUCvA1q5wpz4zFdj3K7D9ayDhgnbMxRto+TzQ9GnApQwvJVEpMaugQqZfi7L0QDTeX3YE8amZcLCzwctdarAWhfIXuQX44ykgJU4LAY/OBKp3LPjZsrUFglpqW89PtNE30jxzdBmQclVrSpLNvSJQp68WWvybat93N+mJ2vf8M1ULPkK+X5p3Gj9lmiOPiCwcgwoVi9jENLyzOBzrjmq1KGGVPVUtSmhFNslRPs03O6YC68YBeh1QsR7Q/zegbFDRT5XUmFRrr229PteGCkt/lmN/aR1zd36vbZ7+WmiR5iG/RrdqblKvAzunAzunabeFVwDQ5mWgwWDOJEtkrUElKSkJp06dMuxHRkbiwIED8Pb2RmAge8+bi6UHLqpalBspWi3Ki51C8GyH6nDgiB66XUYysGy01lwjZFSPrHvj4FJ858rOQRsNJFvWl8Dpv7WalmMrtWYcCUmyla2izXcidv0EZCRqt72ra3Og1Htcey4ist7OtJs2bULHjndW9Q4dOhSzZ88u1s44VPwysnT44K8jmLvzvNqv46fVotSqxGtB+ZDJ2xYMBmKPanOO9JgMNB159/4oJdFpV4Y6S0dcmT02MyXv/eVrawFFwkth+7QQUaEU5vPbZEb9FAWDinE7zD43Zx/2nruuPmdGdwrB6E7BrEWh/J1YCywaCaTFA+4VgMd+1vqWGLNm58QabeSQBBYZZlyj5737rxBRsbHYUT9kGiScPPfbXsQmpsPT2R5fD2yIjjXLG7tYZIpkgjSZHG3jx9oCff7NgMd/ATwrGbdcjm5aPxXZiMikMahQoczbdR7jloYjM1uPGhXcMX1IE1TxdeNZpDtJ7cniZ4HjK7X9JiO05p6SmpSNiCwSgwoVqT9Kz7CK+Oyx+nB34n8hykfsMWDBICDulDbF/INfAA0H81QRUaHxU4YK3R/ltW418XyH6mo2YaI7HFkCLHkeyEzWhgP3/xWo3IgnioiKhEGF7on9UajAdNnAhgnA9q+0/artgEdnAW6+PIlEVGQMKnRX7I9CBZZyDfhzGHBmk7YvM7l2/qB41tshIqvGdxH6z/4oPepUxOePsz8K3UX0AWDBEG1VY1lRuM9UIOwRni4iKhYMKpQH+6NQoRyYByx/GchKA8pWBQbMASrU4UkkomLDoEL59kfxcLbHlAEN0TGU86NQPrIygLXvALuma/sh3YB+P3JVYSIqdgwqpMzfdR7v3ZwfJaS8O6Y/2QRVOT8K5chIAS4dBqL3a9v5HcCNc9p97d8E2o/lrK5EVCIYVKyc9EcZ/9cRzGF/FMq9Js7lI0D0Pq3/iQSTKxHaSse5OXkCD/8AhPbiuSOiEsOgYsXYH4VUE44sEphTUyKb7Ouy7jw5skaPX8NbW0AzwKUsTyIRlSgGFSu177zWH+VyAvujWI3sTODKsbyhRGpOsjPufKyrD+DXKFcwaQB4VCq9lY6JiG5iULHS/ijjlh5BRraO/VEsefK1qyfyhhLpYyKjc27nXCZvTYlsXv4MJURkEhhUrAj7o1iwzFTgwm7g7Hbg7Datf0lmyp2Pk34llernDSVlqzCUEJHJYlCxEuyPYoHBJGqXFkrObddCyu1NOA5ud4YS72ocnUNEZoVBxQqwP4qFDA++cDOYSK3JxT13BhP3ikCVNtoW2ALwrQHY2hmrxERExYJBxcKxP4qZyki+VWMi28W9gC4z72Okc2tOMKnSVqstYWdXIrIwDCoWSqfTY9KqCPy4NVLtc70ecwgmO3MFk335BBO/XMGkDYMJEVkFBhULlJaZjTG/H8DKw5fU/piuNTC6UzBs+Ne26UhPyhtMpPPr7XOXeFbWakqqtNaCiaylw2tIRFaGQcXCXEvOwNO/7FHr9jja2eKzx+qhT4PKxi4WyXBh6fB6ci0QuUUbLnx7MPEKuFVbEtSao3GIiBhULMu5uGQ8NWs3Iq8mw9PZXq3X06Kaj7GLZb2S44DTG4ATa7Svqdfz3u8VmLcpp2yQsUpKRGSyWKNiIfafv46RP+9BXHIGKpdxwc/DmyK4vIexi2Vd9HptUrWTa4ATa7WRObnXx5GJ1YI7A9U7M5gQERUQg4oFWHPkEl6avx9pmTqEVfbEzKeaoryHs7GLZR3SE4Ezm7VwcnIdkBiT9/4KYUBIVyCkO+DfFLDjrxwRUWHwXdPMzdoeiQnLj6o/5juFlsc3AxvCzYmXtURdPaX1NZFwInOa5B6d4+AKVOtwM5x006aiJyKiIuMnmhkPP564MgIztmnDjwc1D8T43nVgb2dr7KJZnqx0bWSOCidrgWtn8t4vo3FqdNfCSVAbwIG1WURExYVBxUyHH7+y4ABWhWvDj9/sEYpn21fj8OPiFH/xVjCRpp3M5Fv32ToAQa1uhpNugE8whw0TEZUQBhUzHH488ufd2Hf+BocfF6esDK3zq/QzkXByOfzO6emlxkTCSdX2gLNnsb48ERHlj0HFjJy9KsOPd+FsXIoafvzjk03QnMOPiz6vyaVDWm2JzGtyfsdtqw3bAP5NtE6wNboBFeux1oSIyAgYVMyETOAmE7lJjYp/WRfMHsbhx4UivY2vHNNCiWxntwJp8Xkf4+pzsyNsd20YsZtvsV5DIiIqPAYVM7A6XBt+nJ6lQ93KXpjxVBMOPy5IMLl+Foi8WWMSuRVIjs37GCdPbQbYqu20rXxtwJadkYmITAmDiombuS0SH67Qhh93Di2PKRx+fHcJ0Vogyak1iT+f9357FyCwxc1g0h6oVJ/zmhARmTgGFROVLcOPV0Rg5nZt+PHgFoH44CEOP84j5ZrWhJPTzyTuZN77be21SdZyakzktr1TKV5FIiK6XwwqJjr8+OX5B7D6iDb8eGzPUPyvHYcfIy1B6/Sqakw2a9PV52ED+DW4FUwCWwKObka5hkREVDwYVExMXFI6Rv6yB/tvDj/+v8fr46H6frBa188BEcuAiOXa6sP67Lz3S7+SnGAic5u4lDVWSYmIqAQwqJiQyJvDj8/FpcDLxUENP25W1RtWJ+40cHSptsUcuHMW2JxgIpt7eWOVkoiISgGDigkNP5aJ3K6nZN4cftwMweXdYTVij2nBRGpPck+2ZmOrjcyp3UebbK1MoDFLSUREpYxBxQSsOhyDlxccUMOP6/l7YcbQpijnYeGdPmUYkwSSnJqTqyfydoKV2pJavYHQBwH3csYsKRERGRGDipHJooIf3Rx+3KWWNvzY1dFCL4v8kNH7gKPLtHByXRvRpNg5AtU6ajUnNXsCrlbY5EVERHew0E9E05OZrUN8amaebUPEZfz2rzbXx5AWQfigdx3Y2drAouh0WifYnGad+Khb99k7A8FdbjXrOHsZs6RERGSCGFQKOWw44WbISEjLFTpS5GtWnhAi9+c8VraUjNtGq+TyVs9QPGNJw49lHZ1z/9wcrfMXkBhz6z4HN23tHGnWkZWHnayoHw4RERUag0o+/jl1VTXJ3F4DIn1I7peHs70a0ePp7ABvN0cMaRmE7nUqwuxlZ2qTr0nNybEVQPKVvFPV1+ih1ZzIGjoOLsYsKRERmREGlXzEJWdgw7Hb1oW5SSo9JGRI2Mi9eea5rYWR2zcPZwfLadqRydekM2zMIa3fycm1QOr1W/c7l9E6wtburS30xxlhiYioCBhU8tEgoAwm96t7ZwhxdYC7oz1sLSVsFLQDrDTdyCywEkou3dxkwb/bufoCtR7UmnVk1I6dgzFKTEREFoRBJR8B3q4Y0MwK5+uQviUy2VpOGFHB5DCQcjX/x3v6A5XqARXCgKptgcBWXOSPiIiKFYOKtcpIAWIjgEsHb9WWXD4CZKXe+ViZdM23JlCxrhZM5GvFehxCTEREJY5BxRokx92qJZFQIptMsKbPp3Owg6tWQ6LCyM1gIuvpsAMsEREZAYOKpcjOAm6cA+JO5d2ungISo/P/HulTompIbtaSVKoPeFcDbO1Ku/RERET5YlAxu46tl24LI6e1rzLLqy7r7t8rASSnyUY2CSjuFbRhTERERCaKQcUUpcXnDSG5Q0lG0t2/z94F8Kl+cwsBfIK12+VCAWfP0vwJiIiIigWDijFrRyR8XDl+Z+1Icv5zuBg6tpYJ0kKIb8jNUCKBJBjw8ANsbUvzpyAiIipRDCqlKfUGELkZOLkOOLXh7n1HhDTL5ASQ3FvZKoC9Y2mWmoiIyGgYVEp6QT4ZaXNqvRZMonYC+uy8i/JJs4yhduRmU413dTbVEBERMaiUgJRrwOm/b4WT25txpO9ISFdtzZug1hz2S0REdA+sUSmO2Vyj999szlkPXNwrHVDyrhZcrT0Q3EULJ9J0Q0RERAXCoFIUSbFabcmpdVrtSe7F+ET5OlookZqTgBbsU0JERFREDCoFkZ0JXNit1ZhIzYn0O8nNyQuo3uFmrUkXwNOvqNeDiIiIcmFQuZv4izf7mawHzmwG0uPz3l+pwa1g4t+Ui/ERERGVAAaV/Pz7PbD6zbzHXLy15hwJJtU7Ae7lS+J6EBERUS4MKvmp3EibWK1yYyBYRuh0AfwacA0cIiKiUsagkh8JKK+fBly9S/t6EBERUS6cbz0/snowQwoREZHRMagQERGRyWJQISIiIpPFoEJEREQmi0GFiIiITBaDChEREZkskwgq3377LapUqQJnZ2c0b94cu3btMnaRiIiIyAQYPagsWLAAY8aMwfvvv499+/ahfv366N69O2JjY41dNCIiIrL2oPLFF1/g6aefxrBhw1C7dm18//33cHV1xcyZM41dNCIiIrLmoJKRkYG9e/eiS5cutwpka6v2d+zYccfj09PTkZCQkGcjIiIiy2XUoHL16lVkZ2ejQoUKeY7L/qVLl+54/KRJk+Dl5WXYAgICSrG0REREZHVNP4Xx1ltvIT4+3rBFRUUZu0hERERkqYsS+vr6ws7ODpcvX85zXPYrVqx4x+OdnJzURkRERNbBqDUqjo6OaNy4MTZs2GA4ptPp1H7Lli2NWTQiIiKy9hoVIUOThw4diiZNmqBZs2b46quvkJycrEYB/Re9Xq++slMtERGR+cj53M75HDfpoNK/f39cuXIF48aNUx1oGzRogNWrV9/RwTY/iYmJ6is71RIREZkf+RyXwTH3YqMvSJwxUdJMFB0dDQ8PD9jY2BR72pMAJB12PT09i/W5idfD3PH3w7TwepgWXo//JtFDQoqfn5+alsSka1Tuh/xw/v7+JfoaElIYVEwHr4dp4fUwLbwepoXX497+qybFLIcnExERkXVhUCEiIiKTxaByFzJfiyyUyHlbTAOvh2nh9TAtvB6mhdejeJl1Z1oiIiKybKxRISIiIpPFoEJEREQmi0GFiIiITBaDChEREZksBpV8fPvtt6hSpQqcnZ3RvHlz7Nq1q/SvjIWZNGkSmjZtqmYRLl++PPr27Yvjx4/neUxaWhpGjRoFHx8fuLu745FHHrljZe3z58/jgQcegKurq3qe119/HVlZWXkes2nTJjRq1Ej1vA8ODsbs2bNL5Wc0Z5MnT1azO7/88suGY7wepe/ixYsYPHiw+h1wcXFB3bp1sWfPHsP9MvZBlhupVKmSur9Lly44efJknue4du0aBg0apCYbK1OmDEaMGIGkpKQ8jzl06BDatm2r3uNkBu5PP/201H5Gc5GdnY333nsPVatWVee6evXq+PDDD/OsTcPrUUpk1A/dMn/+fL2jo6N+5syZ+iNHjuiffvppfZkyZfSXL1/maboP3bt318+aNUsfHh6uP3DggL5Xr176wMBAfVJSkuExzz77rD4gIEC/YcMG/Z49e/QtWrTQt2rVynB/VlaWPiwsTN+lSxf9/v379StXrtT7+vrq33rrLcNjzpw5o3d1ddWPGTNGf/ToUf0333yjt7Oz069evZrX7y527dqlr1Klir5evXr6l156idfDSK5du6YPCgrSP/XUU/qdO3eq/8tr1qzRnzp1yvCYyZMn6728vPRLlizRHzx4UN+7d2991apV9ampqYbH9OjRQ1+/fn39v//+q9+6das+ODhYP3DgQMP98fHx+goVKugHDRqkfh/nzZund3Fx0f/www+l/jObsokTJ+p9fHz0y5cv10dGRur/+OMPvbu7u/7rr782PIbXo3QwqNymWbNm+lGjRhn2s7Oz9X5+fvpJkyaV0iWxDrGxsfJniX7z5s1q/8aNG3oHBwf1ZpAjIiJCPWbHjh1qX4KJra2t/tKlS4bHTJs2Te/p6alPT09X+2+88Ya+Tp06eV6rf//+KijRnRITE/UhISH6devW6du3b28IKrwepe/NN9/Ut2nT5q7363Q6fcWKFfWfffaZ4ZhcJycnJxU2hIRz+Z3ZvXu34TGrVq3S29jY6C9evKj2v/vuO33ZsmUNvzM5r12zZs0S+snM0wMPPKAfPnx4nmP9+vVTAU/wepQeNv3kkpGRgb1796rq1NzrCcn+jh07SquSyyrEx8err97e3uqrnPfMzMw85z40NBSBgYGGcy9fpSo898ra3bt3VwuAHTlyxPCY3M+R8xhev/xJU5s0pd1+zng9St+yZcvQpEkTPPbYY6pZs2HDhvjxxx8N90dGRqoV5nNfK1krRZqnc/+OSHOPPE8Oeby8j+3cudPwmHbt2sHR0THP74g0xV6/fr2UflrT16pVK2zYsAEnTpxQ+wcPHsS2bdvQs2dPtc/rUXrMelHC4nb16lXVLpn7g1DI/rFjx4xWLksjq15LX4jWrVsjLCxMHZM3YHnjlDfZ28+93JfzmPyuTc5993qMhJnU1FTV1kya+fPnY9++fdi9e/cdp4TXo/SdOXMG06ZNw5gxY/D222+r6/Liiy+q34uhQ4ca/o/n9/879/9/CTm52dvbqz8Icj9G+l3c/hw595UtW7ZEf05zMXbsWPW+IX8w2dnZqc+GiRMnqv4/gtej9DCokFH+ig8PD1d/nZBxREVF4aWXXsK6detUh0oyjQAvNSEff/yx2pcaFfk9+f7771VQodL1+++/Y86cOZg7dy7q1KmDAwcOqD+w/Pz8eD1KGZt+cvH19VXJ+faRJrJfsWLF0r42FumFF17A8uXLsXHjRvj7+xuOy/mVprcbN27c9dzL1/yuTc5993qMjIBgbUrepp3Y2Fg1Okr+4pZt8+bNmDJlirotf2HzepQuGclTu3btPMdq1aqlRrrl/j9+r/cn+SrXNTcZFScjgQrze0RQIwqlVmXAgAGqyXnIkCF45ZVX1AhGXo/SxaCSi1SxNm7cWLVL5v4rR/ZbtmxZypfGskjHbQkpixcvxt9//31H1bOcdwcHhzznXtrM5U0659zL18OHD+d5I5YaAQkhOW/w8pjcz5HzGF6/vDp37qzOpfyVmLPJX/NSrZ1zm9ejdElT6O1D9qV/RFBQkLotvzMSJHL//5amCel7kvt3RMK+BNEc8vsm72PSlyXnMVu2bFF9wnL/jtSsWZPNPrmkpKSovj25yR+yci55PUpZKXbcNZvhydKLfvbs2aoH/TPPPKOGJ+ceaUKF99xzz6lhlZs2bdLHxMQYtpSUlDzDk2XI8t9//62GJ7ds2VJttw9P7tatmxriLEOOy5Url+/w5Ndff12NGvr22285PLmAco/64fUwzjBxe3t7NSz25MmT+jlz5qj/y7/99lue4bDyfrR06VL9oUOH9H369Ml3eHLDhg3VEOdt27apUV25hyfLSCEZnjxkyBA1PFne8+R1ODw5r6FDh+orV65sGJ68aNEiNR2CjCzk9ShdDCr5kLk35ANT5lOR4coyHwHd5380IN9N5lbJIW+2zz//vBo6KW+cDz/8sAozuZ09e1bfs2dPNe+DvGm8+uqr+szMzDyP2bhxo75Bgwbq+lWrVi3Pa1DBgwqvR+n766+/VBiXP5ZCQ0P106dPz3O/DIl97733VNCQx3Tu3Fl//PjxPI+Ji4tTwUTm/JCh+8OGDVPD0HOTOVhkKLQ8h3wYSwCivBISEtTvg3wWODs7q/eSd955J8+wbl6P0mEj/5R2LQ4RERFRQbCPChEREZksBhUiIiIyWQwqREREZLIYVIiIiMhkMagQERGRyWJQISIiIpPFoEJEREQmi0GFiIiITBaDCpGVmz17NsqUKXPfz9OhQwe1umxpeeqpp9C3b1+YgrNnz8LGxkatk0RExYtBhcjK9e/fXy1+R+YXkIisgb2xC0BExuXi4qI2IiJTxBoVIguzfPly1ZSTnZ2t9qU5Qpolxo4da3jMyJEjMXjw4Hybfj744AM0aNAAv/76K6pUqQIvLy8MGDAAiYmJhsckJyfjySefhLu7OypVqoT/+7//u6Mc169fV48pW7YsXF1d0bNnT5w8eVLdJ0uMlStXDn/++afh8fKa8lw5tm3bBicnJ6SkpBTo59bpdJg0aRKqVq2qglf9+vXzPP+mTZvUediwYQOaNGmiytSqVSscP348z/N89NFHKF++PDw8PNR5kvMmZcs5Nz///DOWLl2qnks2ed4cZ86cQceOHdVzy+vv2LGjQGUnortjUCGyMG3btlWhYv/+/Wp/8+bN8PX1zfOBKsekT8ndnD59GkuWLFGhRzZ5/OTJkw33v/766+qYfGCvXbtWPfe+ffvuaCLZs2cPli1bpj6wJZz06tULmZmZ6gO+Xbt2hjJJqImIiEBqaiqOHTtmKGPTpk3Vh35BSEj55Zdf8P333+PIkSN45ZVXVBiT58ntnXfeUcFKymZvb4/hw4cb7pszZw4mTpyITz75BHv37kVgYCCmTZtmuP+1117D448/jh49eiAmJkZtEnZyP7c8RsJhjRo1MHDgQGRlZRWo/ER0F6W0SjMRlaJGjRrpP/vsM3W7b9+++okTJ+odHR31iYmJ+gsXLsiK6foTJ06o+2fNmqX38vIyfO/777+vd3V1Vcvc53j99df1zZs3V7flOeS5fv/9d8P9cXFxehcXF/1LL72k9uW55TW2b99ueMzVq1fVY3K+b8qUKfo6deqo20uWLFHP36dPH/20adPUsS5duujffvvtu/6MQ4cOVY8XaWlpqsz//PNPnseMGDFCP3DgQHV748aNqkzr16833L9ixQp1LDU1Ve1LGUaNGpXnOVq3bq2vX79+vq+bIzIyUj3PTz/9ZDh25MgRdSwiIuKuPwMR/TfWqBBZoPbt26vaCqnF2Lp1K/r164datWqp5hSpYfDz80NISMhdv1+afKTpI4c0ycTGxhpqWzIyMtC8eXPD/d7e3qhZs6ZhX2pHpLYi92N8fHzUY+S+nDIePXoUV65cMdTwyCblllqXf/755561PrmdOnVKNRF17dpVNUflbFLDIuXNrV69enl+LpHzs0kzULNmzfI8/vb9e7nXcxNR0bAzLZEFkg/4mTNn4uDBg3BwcEBoaKghBEgzi4SEe5HvyU2aaqQPSHGqW7euCjgSUmSTJpeKFSuqZpfdu3ersJK7WeVekpKS1NcVK1agcuXKee6Tfi53+9nk5xLF9bOV5HMTWSvWqBBZcD+VL7/80hBKcoKKbAWtqchP9erV1Qfyzp07Dcck/OQe4iy1N9I3I/dj4uLiVI1F7dq1DR/kUk7p5yJ9Stq0aaNqJNLT0/HDDz+oDq9ubm4FKpM8pwSS8+fPIzg4OM8WEBBQ4J9NanwkJOV2+76jo6OhozIRlTzWqBBZIBlpIx/60jl06tSp6ph0XpWOoFJT8V81KvciTSojRoxQHWqlOUdGyEgnUlvbW3/3SLNSnz598PTTT6vQIc1IMnpGajvkeA4JTK+++qoKJfK8OeWUcsvzF5Q8v3RilQ60UoMhoSc+Ph7bt2+Hp6cnhg4dWqDnGT16tCqzlEdqcxYsWIBDhw6hWrVqeZrF1qxZo0KX/PwyKoqISg5rVIgslIQR+cs/p/ZEmlmk5kGaV3L3JymKzz77TNWGPPTQQ+jSpYsKBo0bN87zmFmzZqljDz74IFq2bKn6y6xcuTJP88jtZRRy+/ZjBfHhhx/ivffeU6N/pEZHRuZIU5AMVy6oQYMG4a233lKhp1GjRoiMjFSjl5ydnQ2PkSAj50/CjAyxljBERCXHRnrUluDzExGZNemgK+FO5pUhotLHph8ioptk5JDMw9K9e3fY2dlh3rx5WL9+PdatW8dzRGQkrFEhIrpJJpyT5iyZLC8tLU018bz77rtqeDcRGQeDChEREZksdqYlIiIik8WgQkRERCaLQYWIiIhMFoMKERERmSwGFSIiIjJZDCpERERkshhUiIiIyGQxqBARERFM1f8DcPveoIkRaO0AAAAASUVORK5CYII=", "text/plain": [ - "
" + "
" ] }, - "metadata": { - "needs_background": "light" - }, + "metadata": {}, "output_type": "display_data" } ], @@ -470,7 +468,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.10.20" } }, "nbformat": 4,