diff --git a/melange/csmc.py b/melange/csmc.py index a018bc7..5614f55 100644 --- a/melange/csmc.py +++ b/melange/csmc.py @@ -4,18 +4,21 @@ from melange.smc_objects import StaticULA from melange.gaussians import * import jax.numpy as jnp +from jax import ops from functools import partial from jax import random -from jax import grad, vmap, jit +from jax import grad, vmap, jit, device_put from jax.lax import map, scan, cond from jax.scipy.special import logsumexp from jax.config import config; config.update("jax_enable_x64", True) +import numpy as np +import tqdm import logging # Instantiate logger logging.basicConfig(level = logging.NOTSET) _logger = logging.getLogger("csmc") -_logger.setLevel(logging.INFO) +_logger.setLevel(logging.DEBUG) """twisting utilities""" def twisted_normal_params(mu, cov, A, b, get_log_twist_constant=True): #tested @@ -119,7 +122,15 @@ def do_param_twist(base_mu, base_cov, As, bs, get_log_normalizer): logZ = cond(get_log_normalizer, lambda x: Normal_logZ(*x) - base_logZ, lambda x: 0., (twisted_mu, twisted_cov)) # compute the logZ of the kernel if specified return twisted_mu, twisted_cov, logZ + log_twist_constant -def gaussians_twist(Xp, potential, dt, potential_params, A_fn, b_fn, A_params, b_params, get_log_normalizer): +def gaussians_twist(Xp, + potential, + dt, + potential_params, + A_fn, + b_fn, + A_params, + b_params, + get_log_normalizer): """ conduct a gaussian twist to return a twisted mean, covariance vector, normalizing constant, As, and bs @@ -161,6 +172,44 @@ def gaussians_twist(Xp, potential, dt, potential_params, A_fn, b_fn, A_params, b twisted_mu, twisted_cov, logZ = do_param_twist(start_mu, start_cov, As, bs, get_log_normalizer) return twisted_mu, twisted_cov, logZ, (As, bs) +def precomputed_gaussian_twist(Xp, + twisted_mu, + twisted_cov, + A_fn, + b_fn, + A_params, + b_params): + """ + conduct a gaussian twist with recomputed twisted mu and covariance vectors + NOTE : this is primarily a utility for ADP + + arguments + Xp : jnp.array(Dx) + previous x values + twisted_mu : jnp.array(Dx) + twisted mean + twisted_cov : jnp.array(Dx) + twisted covariance vector + A_fn : function + A twisting function + b_fn : function + b twisting function + A_params : jnp.array(V) + parameters to twist; each is of dim V + b_params : jnp.array(W) + parameters to twist; each is of dim W + returns + logZ : float + forward kernel log normalizer + """ + A, b = A_fn(Xp, A_params), b_fn(Xp, b_params) + _, _, logZ = do_param_twist(twisted_mu, + twisted_cov, + A_params[jnp.newaxis, ...], + b_params[jnp.newaxis, ...], + True) + return logZ + def get_twisted_gmm(mix_weights, mus, covs, A, b): """ twist a gaussian mixture model with matrix A and vector b @@ -201,13 +250,94 @@ def get_twisted_gmm(mix_weights, mus, covs, A, b): return full_log_normalizer, log_normalized_twisted_mixtures, (twisted_mus, twisted_covs) -"""gaussian twisting potentials""" -log_psi_twist_A = lambda x, As: -(x*(As.sum(axis=0))).dot(x) -log_psi_twist_b = lambda x, bs: -x.dot(bs.sum(axis=0)) -log_psi_twist = lambda x, As, bs: log_psi_twist_A(x, As) + log_psi_twist_b(x, bs) +""" +gaussian twisting potentials + +log_psi_twist_A and log_psi_twist_b are +""" +def log_psi_twist_A(x, As): + """ + compute the log twisted potential of the A values + arguments + x : jnp.array(Dx) + position + As : jnp.array(Q, Dx) + Q covariance twisting matrices + + returns + out : float + negative square mahalonobis distance + """ + return -(x*(As.sum(axis=0))).dot(x) + +def log_psi_twist_b(x, bs): + """ + compute the log twisted potential of the b values + + arguments + x : jnp.array(Dx) + position + bs : jnp.array(Q, Dx) + Q mean twisting vector + + returns + out : float + twisted b potential + """ + return -x.dot(bs.sum(axis=0)) + +def log_psi_twist(x, As, bs): + """ + compute the full log twist + + arguments + x : jnp.array(Dx) + position + As : jnp.array(Q, Dx) + Q mean twisting matrices (diagonal) + bs : jnp.array(Q, Dx) + Q mean twisting vector + + return + out : float + full twist + """ + return log_psi_twist_A(x, As) + log_psi_twist_b(x, bs) + +""" +define a vlog_psi_twist to vectorize all of the log_psi_twist inputs + + arguments + xs : jnp.array(N, Dx) + vectorized positions + N_As : jnp.array(N, Q, Dx) + vectorized covariance twists + N_bs : jnp.array(N, Q, Dx) + vectorized twisting vectors + return + out : jnp.array(N) +""" vlog_psi_twist = vmap(log_psi_twist, in_axes=(0,0,0)) -vlog_psi_twist_mtnr = vmap(log_psi_twist, in_axes=(0,None, None)) +#this util is not used directly +vlog_psi_twist_util = vmap(log_psi_twist, in_axes=(0,None, None)) + +def vlog_psi_twist_single(xs, A, b): + """ + do a twist on N different positions with a single A and b vector + + arguments + xs : jnp.array(N, Dx) + vectorized positions + A : jnp.array(Dx) + vectorized covariance twists + b : jnp.array(Dx) + vectorized twisting vectors + return + out : jnp.array(N) + + """ + return vlog_psi_twist_util(xs, A[jnp.newaxis, ...], b[jnp.newaxis, ...]) ######################################### #####CSMC Classes######################## @@ -227,25 +357,47 @@ def __init__(self, *args, # N, potential, forward_potential, backward_potential A_fn, b_fn, - A_params_len, - b_params_len, + params_len, Dx, T, - max_twists = 10, + max_twists, + prop_params, + init_params, + model_params, + y, **kwargs): _logger.debug(f"instantiating super") super().__init__(*args, **kwargs) #init the og args + #define the twisting functions + _logger.debug(f"define A, b functions") + self.A_fn = A_fn + self.b_fn = b_fn + + _logger.debug(f"getting the parameter lengths") + self.params_len = params_len + _logger.debug(f"getting x dimension") self.Dx = Dx _logger.debug(f"getting the terminal time") self.T = T - #define the twisting functions - _logger.debug(f"define A, b functions") - self.A_fn = A_fn - self.b_fn = b_fn + _logger.debug(f"getting the maximum number of twists") + self.max_twists = max_twists + + _logger.debug(f"getting prop_params") + self.prop_params = prop_params + + _logger.debug(f"getting init_params") + self.init_params = init_params + + _logger.debug(f"getting model_params") + self.model_params = model_params + + _logger.debug(f"getting y") + self.y = y + #vmap the A, b functions; this is a special utility to compute As, bs of a single twist for N particle positions self.vA_fn = vmap(A_fn, in_axes=(0,None)) @@ -254,45 +406,47 @@ def __init__(self, #define twisting iterations _logger.debug(f"define twisting parameters and parameter caches") self.twisting_iteration = 1 - self.max_twists = max_twists - self.A_params_cache = jnp.zeros((max_twists, A_params_len)) - self.b_params_cache = jnp.zeros((max_twists, A_params_len)) + self.A_params_cache = jnp.zeros((max_twists, T, params_len)) + self.b_params_cache = jnp.zeros((max_twists, T, params_len)) self.A0, self.b0 = jnp.zeros(self.Dx), jnp.zeros(self.Dx) - #assert parameter lengths are the same - assert A_params_len == b_params_len, f"the A and b parameter lengths must be the same" - self.A_params_len = A_params_len - self.b_params_len = b_params_len - #define the psi twisting function _logger.debug(f"define psi twisting functions") self.vlog_psi_twist = vmap(log_psi_twist, in_axes=(0,0,0)) + self.precomputed_twist_fn = precomputed_gaussian_twist self.twist_fn = gaussians_twist self.vtwist_fn = vmap(self.twist_fn, in_axes=(0, None, None, None, None, None, None, None, None)) #args (Xp, potential, dt, potential_params, A_fn, b_fn, A_params, b_params, get_log_normalizer) + self.vprecomputed_twist_fn = vmap(self.precomputed_twist_fn, in_axes=(0, 0, 0, None, None, None, None)) # (Xp, twisted_mu,twisted_cov, A_fn, b_fn, A_params, b_params) + + #make a cache to store the As, bs, logK_twists, twisted_mus, twisted_covs for the next iteration and to compute ADP + self.As_cache = None # cache of the A values + self.bs_cache = None # cache of the b values + self.K_logZs_cache = None # cache of the twisted kernel log normalizers + self.twisted_mus_cache = None # cache of the twisted mus + self.twisted_covs_cache = None # cache of the twisted covariances + self.t0_log_normalizer = None # _value_ of the prior normalizer def sim_prop_fn(self): _logger.debug(f"generating simulation propagation function") self.vpropagator = vmap(gaussian_proposal, in_axes=(0,0,0)) - def prop(t, Xp, y, prop_params, model_params, rs): - N, Dx = Xp.shape #extract the number of particles and the dimension of x - folder_rs = random.split(rs, num=N+1) #create a list of random keys + def prop(t, Xp, rs): + folder_rs = random.split(rs, num=self.N+1) #create a list of random keys new_rs, runner_rs = folder_rs[0], folder_rs[1:] #separate the random keys - twisted_mus, twisted_covs, _, _ = self.vtwist_fn(Xp, - self.forward_potential, - prop_params['dt'], - prop_params['forward_potential_params'][t], - self.A_fn, - self.b_fn, - self.A_params_cache[:self.twisting_iteration], - self.b_params_cache[:self.twisting_iteration], - False) # twisted_mu, twisted_cov, logZ, (As, bs) - + twisted_mus, twisted_covs, K_logZs, (As, bs) = self.vtwist_fn(Xp, + self.forward_potential, + self.prop_params['dt'], + self.prop_params['forward_potential_params'][t], + self.A_fn, + self.b_fn, + self.A_params_cache[:self.twisting_iteration,t], + self.b_params_cache[:self.twisting_iteration,t], + True) # twisted_mu, twisted_cov, logZ, (As, bs) Xs = self.vpropagator(runner_rs, twisted_mus, twisted_covs) - return Xs + return Xs, {'twisted_mus': twisted_mus, 'twisted_covs': twisted_covs, 'K_logZs': K_logZs, 'As': As, 'bs': bs} return prop def log_weights_fn(self): @@ -300,7 +454,7 @@ def log_weights_fn(self): #define vmapped twisting functions and logK_integral functions static_ula_logW_fn = super().log_weights_fn() #get the super static ula log_weights_fn - def log_weights(t, Xp, Xc, y, prop_params, model_params): + def log_weights(t, Xp, Xc): """ arguments prop_params : tuple @@ -309,30 +463,16 @@ def log_weights(t, Xp, Xc, y, prop_params, model_params): backward_potential_params : jnp.array(T-1,S) forward_dts : jnp.array(T) backward_dts : jnp.array(T-1) - - NOTE : need to implement a conditional here. specifically, if t == T, then we do _not_ need to compute the logK_Zs. we need only the untwisted incremental weight and the - twisting potential """ # compute untwisted logWs - logWs = static_ula_logW_fn(t, Xp, Xc, y, prop_params, model_params) - - #compute twisting logZ forward kernel and psi_twist - compute_logK_Z = cond(t == self.T, lambda x: False, lambda x: True, None) - _, _, logK_Zs, (As, bs) = self.vtwist_fn(Xp, - self.forward_potential, - prop_params['dt'], - prop_params['forward_potential_params'][t], - self.A_fn, - self.b_fn, - self.A_params_cache[:self.twisting_iteration], - self.b_params_cache[:self.twisting_iteration], - compute_logK_Z) + logWs = static_ula_logW_fn(t, Xp, Xc, self.y, self.prop_params, self.model_params) #compute twisting functions - log_psi_ts = self.vlog_psi_twist(Xc, As, bs) + log_psi_ts = self.vlog_psi_twist(Xc, self.As_cache[t], self.bs_cache[t]) #build modifier - return logWs + logK_Zs - log_psi_ts + K_logZs = cond(t==self.T-1, lambda x: jnp.zeros(self.N), lambda x: self.K_logZs_cache[t+1], None) #compute the t+1 twisted kernel normalizer + return logWs + K_logZs - log_psi_ts return log_weights def initialize_Xs_fn(self): @@ -344,34 +484,32 @@ def initialize_Xs_fn(self): self.vsample_gmm = vmap(sample_gmm, in_axes = (0, None, None, None)) # args: (key, weights, mus, covs) - def init_xs(prop_params, rs, init_params): + def init_xs(rs): #resolve random keys folder_rs = random.split(rs, self.N+1) rs, run_rs = folder_rs[0], folder_rs[1:] # conduct twisted initialization full_log_normalizer, log_normalized_twisted_mixtures, (twisted_mus, twisted_covs) - _, log_normalized_twisted_mixtures, (twisted_mus, twisted_covs) = get_twisted_gmm(init_params['mixture_weights'], - init_params['mus'], - init_params['covs'], + t0_log_normalizer, log_normalized_twisted_mixtures, (twisted_mus, twisted_covs) = get_twisted_gmm(self.init_params['mixture_weights'], + self.init_params['mus'], + self.init_params['covs'], self.A0, self.b0) #compute twisting parameters _, Xs = self.vsample_gmm(run_rs, exp_normalize(log_normalized_twisted_mixtures), twisted_mus, twisted_covs) - return Xs + return t0_log_normalizer, Xs return init_xs def initialize_logW_fn(self, **kwargs): _logger.debug(f"generating SMC log weight initializer function") - def init_logWs(X, init_params, prop_params): - log_normalizer, _, _ = get_twisted_gmm(init_params['mixture_weights'], - init_params['mus'], - init_params['covs'], - self.A0, - self.b0) #compute twisting parameters - + def init_logWs(X): # compute twisted weights - log_psi0s = log_normalizer - vlog_psi_twist_mtnr(X, self.A0[jnp.newaxis, ...], self.b0[jnp.newaxis, ...]) + log_psi0s = ( + self.t0_log_normalizer # use the precomputed t=0 log normalizer; NOTE : this isn't necessary for ADP + + self.K_logZs_cache[1] # use the precomputed t=0 forward kernel normalizer + - vlog_psi_twist_single(X, self.A0, self.b0) # compute the t=0 twisted potential + ) return log_psi0s return init_logWs @@ -383,18 +521,55 @@ def get_ADP_fn(self): """ _logger.debug(f"generating approximate dynamic programming functions") from jax.ops import index, index_add, index_update - import tqdm import numpy as np def sum_square_diffs(Xcs, As, bs, V_bars): """ compute the sum of square differences between the potential (-vlog_psi_twist) and the precomputed V_bars + + arguments + Xcs : jnp.array(N, Dx) + current positions + As : jnp.array(N, Dx) + A twists + bs : jnp.array(N, Dx) + b twists + V_bars : jnp.array(N) + previously-computed potentials + + returns + out : float + sum of square differences + + NOTE : vlog_psi_twist takes As, bs of dim (N, Q, Dx), but the As, bs have dim (N, Dx); + since we are summing over Q=1, we have to expand the dimension of the As, bs at axis 1 """ - return ((-self.vlog_psi_twist(Xcs, As, bs) - V_bars)**2.).sum() + return ( + (-self.vlog_psi_twist(Xcs, + jnp.expand_dims(As, axis=1), + jnp.expand_dims(bs, axis=1) + ) + - V_bars)**2.).sum() def loss(Xps, Xcs, A_params, b_params, V_bars): """ the loss function is the sum_square_diffs explicitly parameterized by A_params, b_params + + arguments + Xps : jnp.array(N, Dx) + previous positions + Xcs : jnp.array(N, Dx) + current positions + A_params : jnp.array(Q) + parameters to pass to the A twisting function + b_params : jnp.array(Q) + parameters to pass to the b twisting function + V_bars : jnp.array(N) + previously-computed potentials + + returns + out : float + sum of square differences """ As, bs = self.vA_fn(Xps, A_params), self.vb_fn(Xps, b_params) return sum_square_diffs(Xcs, As, bs, V_bars) @@ -402,33 +577,76 @@ def loss(Xps, Xcs, A_params, b_params, V_bars): def scipy_loss(A_b_params, Xps, Xcs, V_bars): """ rewrite the loss function s.t. it is amenable to scipy.optimize.minimize library + + arguments + A_b_params : jnp.array(Q + R) + flattened A, b parameter inputs where Q is the dimension of A params and R is dimension of R params + + see above functions for redundant arguments + + returns + out : float + loss """ - A_params, b_params = A_b_params[:self.A_params_len], A_b_params[self.A_params_len:] + A_params, b_params = A_b_params[:self.params_len], A_b_params[self.params_len:] return loss(Xps, Xcs, A_params, b_params, V_bars) def t0_scipy_loss(A0_b0, Xcs, V_bars): """ the scipy loss function for the t=0 twisting iteration is slightly different since we need not compute a tensor of As, bs. the twisting A0, b0 arrays are specified _exclusively_ as x-independent and parameter-independent arrays + + arguments + A0_b0 : jnp.array(2*Dx) + flattened A0, b0 input parameters + + see above functions for redundant arguments + + returns + out : float + loss """ A0, b0 = A0_b0[:self.Dx], A0_b0[self.Dx:] - return ((-vlog_psi_twist_mtnr(Xcs, A0[np.newaxis, ...], b0[np.newaxis, ...]) - V_bars)**2).sum() + return ((-vlog_psi_twist_single(Xcs, A0, b0) - V_bars)**2).sum() + + def get_callback_fn(Xps, Xcs, V_bars, t0=False): + """ + utility to retrieve a callback function to pass to the scipy.optimize.minimize function for the purpose of logging outputs - def get_callback_fn(Xps, Xcs, V_bars): + arguments + see above functions for redundant arguments + + returns + loss_logger : list() + empty list to log the losses + callback : fn + callback function that is supposed to be passed to the scipy.optimize.minimize function + + """ loss_logger = [] - partial_scipy_loss = partial(scipy_loss, Xps=Xps, Xcs=Xcs, V_bars=V_bars) + if t0: + partial_scipy_loss = partial(t0_scipy_loss, Xcs=Xcs, V_bars=V_bars) + else: + partial_scipy_loss = partial(scipy_loss, Xps=Xps, Xcs=Xcs, V_bars=V_bars) + def callback(xk): val=partial_scipy_loss(xk) loss_logger.append(val) return loss_logger, callback + #jax grad wrapper + scipy_loss_grad = grad(scipy_loss, argnums=(0,)) + t0_scipy_loss_grad = grad(t0_scipy_loss, argnums=(0,)) - + #to numpy wrapper + def np_scipy_loss(*args): return float(scipy_loss(*args)) + def np_t0_scipy_loss(*args): return float(t0_scipy_loss(*args)) + def np_scipy_loss_grad(*args): return np.array(scipy_loss_grad(*args), dtype=np.float64) + def np_t0_scipy_loss_grad(*args): return np.array(t0_scipy_loss_grad(*args), dtype=np.float64) def ADP(Xs, logWs, - prop_params, minimizer_params = { 'method': 'L-BFGS-B', 'options': {'disp':True} @@ -438,6 +656,8 @@ def ADP(Xs, """ conduct ADP + ISSUE : there is a discrepancy between the initial value of the loss function (scipy_loss) when computed inside/outside the scipy's minimize function. + WARNING : this function CANNOT be jit'd """ from scipy.optimize import minimize #import the optimizer @@ -447,91 +667,199 @@ def ADP(Xs, assert T == self.T assert N == self.N - out_A_params = np.zeros((T, self.A_params_len)) - out_b_params = np.zeros((T, self.b_params_len)) - loss_trace=[] - - #make a reporter array - loss_reporter_container = np.zeros((T, 2)) #for each time, we report the initial loss and final loss + #initialize arrays as nan + out_A_params = np.empty((T, self.params_len)); out_A_params[:] = np.nan + out_b_params = np.empty((T, self.params_len)); out_b_params[:] = np.nan + loss_trace=[] #create a loss trace if we are verbose logK_int_twist = jnp.zeros(N) #initialize the logK_twist_T+1 _logger.debug(f"iterating backward...") for t in tqdm.tqdm(jnp.arange(1,T)[::-1]): #do ADP backward Vbar_t = -logWs[t] - logK_int_twist #define the Vbar_ts - init_loss = (Vbar_t**2).sum() - loss_reporter_container[t,0] = init_loss - x0 = np.zeros((2, self.A_params_len)).flatten() # TODO : assertion that A, b params have to same size # TODO : do we have to flatten an array as first argument to minimizer + runner_args = (Xs[t-1], Xs[t], Vbar_t) #other arguments passed to minimizer if verbose: - loss_logger, callback_fn = get_callback_fn(Xs[t-1], Xs[t], Vbar_t) - out_containter = minimize(scipy_loss, - x0 = x0, #initial parameters - args = (Xs[t-1], Xs[t], Vbar_t), #other args + #generally, the loss loggers do not have consistent number of minimization steps at each iteration + #so it is just an empty list + loss_logger, callback_fn = get_callback_fn(*runner_args) + else: + callback_fn=None + + out_containter = minimize(np_scipy_loss, #at time t, we use the generic np scipy loss function + x0 = jnp.zeros((2, self.params_len)).flatten(), #initial parameters + args = runner_args, #other args + jac = np_scipy_loss_grad, callback=callback_fn, **minimizer_params) - _logger.debug(f"minimizer message: {out_containter['message']}") + if verbose: - loss_trace.append(loss_logger) + _logger.debug(f"minimizer message: {out_containter['message']}") + loss_trace.append(np.array(loss_logger)) assert out_containter['success'], f"iteration {t} failed with message {out_containter['message']}" #TODO : make this fail safe - reshaped_out_x = out_containter['x'].reshape(2,self.A_params_len) + + #extract the output parameters + reshaped_out_x = out_containter['x'].reshape(2,self.params_len) #check to make sure this is reshaped properly new_A_params, new_b_params = reshaped_out_x[0,:], reshaped_out_x[1,:] - final_loss = out_containter['fun'] + + #record the parameters out_A_params[t,:] = new_A_params out_b_params[t,:] = new_b_params - loss_reporter_container[t,1] = final_loss #compute the logK twisting integrals for the next optimization step - _, _, logK_int_twist, _ = self.vtwist_fn(Xs[t-1], - self.forward_potential, - prop_params['dt'], - prop_params['forward_potential_params'][t], - self.A_fn, - self.b_fn, - jnp.vstack((self.A_params_cache[:self.twisting_iteration], jnp.asarray(new_A_params))), - jnp.vstack((self.b_params_cache[:self.twisting_iteration], jnp.asarray(new_b_params))), - True) + logK_int_twist = self.vprecomputed_twist_fn(Xs[t-1], + self.twisted_mus_cache[t], # CHECK that this index is correct + self.twisted_covs_cache[t], # Same here + self.A_fn, + self.b_fn, + new_A_params, + new_b_params) + # _, _, logK_int_twist, _ = self.vtwist_fn(Xs[t-1], + # self.forward_potential, + # self.prop_params['dt'], + # self.prop_params['forward_potential_params'][t], + # self.A_fn, + # self.b_fn, + # jnp.vstack((self.A_params_cache[:self.twisting_iteration,t], jnp.asarray(new_A_params))), + # jnp.vstack((self.b_params_cache[:self.twisting_iteration,t], jnp.asarray(new_b_params))), + # True) #now to do the 0th time twist _logger.debug(f"conducting t=0 twist...") - Vbar0 = -logWs[0] - logK_int_twist #compute the Vbars - loss_reporter_container[0,0] = (Vbar0**2).sum() - x0 = np.zeros((2, self.Dx)).flatten() + Vbar0 = -(logWs[0] - self.t0_log_normalizer) - logK_int_twist #compute the Vbars # TODO : do we have to flatten an array as first argument to minimizer - out_containter = minimize(t0_scipy_loss, - x0 = x0, - args = (Xs[0], Vbar0), + + runner_args = (None, Xs[0], Vbar0) #other arguments passed to minimizer + if verbose: + loss_logger, callback_fn = get_callback_fn(*runner_args, t0=True) + else: + callback_fn=None + + out_containter = minimize(np_t0_scipy_loss, + x0 = jnp.zeros((2, self.Dx)).flatten(), + args = runner_args[1:], #we cannot pass 'None' + jac = np_t0_scipy_loss_grad, + callback=callback_fn, **minimizer_params) + if verbose: + _logger.debug(f"minimizer message: {out_containter['message']}") + loss_trace.append(np.array(loss_logger)) assert out_containter['success'], f"iteration 0 failed with message {out_containter['message']}" #TODO : make this fail safe reshaped_out_x = out_containter['x'].reshape(2,self.Dx) new_A0, new_b0 = reshaped_out_x[0,:], reshaped_out_x[1,:] - final_loss = out_containter['fun'] - loss_reporter_container[0,1] = final_loss - _logger.debug(f"finished optimizations") + _logger.debug(f"finished optimizations; terminating...") # package the new parameters - output = {'losses': loss_reporter_container, - 'out_A_params': out_A_params, - 'out_b_params': out_b_params, - 'out_A0': new_A0, - 'out_b0': new_b0, - 'loss_trace': loss_trace} + output = { + 'out_A_params': jnp.asarray(out_A_params[1:]), + 'out_b_params': jnp.asarray(out_b_params[1:]), + 'out_A0': jnp.asarray(new_A0), + 'out_b0': jnp.asarray(new_b0), + 'loss_trace': loss_trace[::-1] + } return output - return {'utilities': (sum_square_diffs, loss, scipy_loss, t0_scipy_loss), 'ADP': ADP} + out_fns = [sum_square_diffs, loss, scipy_loss, t0_scipy_loss, get_callback_fn, scipy_loss_grad, t0_scipy_loss_grad, np_scipy_loss, np_t0_scipy_loss, np_scipy_loss_grad, np_t0_scipy_loss_grad, ADP] + output_fn_dict = {fn.__name__: fn for fn in out_fns} + return output_fn_dict - def update_cSMC(self, new_A_params, new_b_params, new_A0, new_b0): + def update_cSMC(self, + new_A_params, + new_b_params, + new_A0, + new_b0): """ update the instance with the new twisting parameters and initial parameters """ _logger.debug(f"updating cSMC with new parameters") #make sure that the input parameters are compatible with those of the class - assert new_A_params.shape == (self.T, self.A_params_len) - assert new_b_params.shape == (self.T, self.b_params_len) + assert new_A_params.shape == (self.T-1, self.params_len) + assert new_b_params.shape == (self.T-1, self.params_len) assert new_A0.shape == (self.Dx,) assert new_b0.shape == (self.Dx,) - # + #catalogue the parameters + catalogue_index = self.twisting_iteration-1 + _logger.debug(f"cacheing parameters at index {catalogue_index}") + self.A_params_cache = ops.index_update(self.A_params_cache, ops.index[catalogue_index,1:,:], new_A_params) + self.b_params_cache = ops.index_update(self.b_params_cache, ops.index[catalogue_index,1:,:], new_b_params) + self.A0 = self.A0 + new_A0 + self.b0 = self.b0 + new_b0 + + #then update the twisting iteration + self.twisting_iteration += 1 + + +""" +implementation utilities; + +these utilities are used to perform csmc twists, run SIS, etc +""" +def twisted_smc(smc_object, rs, aggregate_works=False): + """ + generate trajectories of length T, N + """ + + #initialize + init_X_fn = smc_object.initialize_Xs_fn() #get the initialize xs function + rs, init_rs = random.split(rs) #split the random vars + t0_log_normalizer, X0 = init_X_fn(init_rs) # get the log normalizing function at t0 and the initial positions + + #make propagation function + prop_fn = jit(smc_object.sim_prop_fn()) #we jit this function because it can become expensive + + #define a cache + As_cache = np.zeros((smc_object.T, smc_object.N, smc_object.twisting_iteration, smc_object.Dx)) # create a cache of A variance parameters + bs_cache = np.zeros((smc_object.T, smc_object.N, smc_object.twisting_iteration, smc_object.Dx)) # createa a cache of b twisting parameters + K_logZs_cache = np.zeros((smc_object.T, smc_object.N)) # create a cache of log normalizing parameters + twisted_mus_cache = np.zeros((smc_object.T, smc_object.N, smc_object.Dx)) # create a cache of twisted mus + twisted_covs_cache = np.zeros((smc_object.T, smc_object.N, smc_object.Dx)) #create a cache of twisted covariance vectors + all_Xs = np.zeros((smc_object.T, smc_object.N, smc_object.Dx)) # create a cache of all positions + all_Xs[0] = np.asarray(X0) #update the t=0 particles + + Xp = X0 # the previous positions are set as the X0s before SMC + for t in tqdm.trange(1,smc_object.T): + rs, run_rs = random.split(rs) #split the random vars + X, out_dict = prop_fn(t, Xp, run_rs) #do a propagation + + #update the positions + As_cache[t] = np.asarray(out_dict['As']) # update the As cache + bs_cache[t] = np.asarray(out_dict['bs']) # and the bs cache + K_logZs_cache[t] = np.asarray(out_dict['K_logZs']) # update the log normalizers + twisted_mus_cache[t] = np.asarray(out_dict['twisted_mus']) # update the twisted mus cache + twisted_covs_cache[t] = np.asarray(out_dict['twisted_covs']) # update the twisted covariance cache + all_Xs[t] = np.asarray(X) # update positions cache + Xp = X + + #update the smc object cache + smc_object.As_cache = jnp.asarray(As_cache) + smc_object.bs_cache = jnp.asarray(bs_cache) + smc_object.K_logZs_cache = jnp.asarray(K_logZs_cache) + smc_object.twisted_mus_cache = jnp.asarray(twisted_mus_cache) + smc_object.twisted_covs_cache = jnp.asarray(twisted_covs_cache) + smc_object.t0_log_normalizer = t0_log_normalizer + all_Xs = jnp.asarray(all_Xs) + + #get weight_functions + init_logw_fn = smc_object.initialize_logW_fn() + logw_fn = jit(smc_object.log_weights_fn()) + + #compute initial weight + init_logWs = init_logw_fn(all_Xs[0]) + + #compute the other weights + logWs = np.zeros((smc_object.T, smc_object.N)) + logWs[0,:] = np.asarray(init_logWs) + + for t in tqdm.trange(1,smc_object.T): + Xp, Xc = all_Xs[t-1], all_Xs[t] + logW = logw_fn(t, Xp, Xc) + logWs[t] = np.asarray(logW) + + logWs = jnp.asarray(logWs) + + + return all_Xs, cond(aggregate_works, lambda x: jnp.cumsum(x, axis=0), lambda x: x, logWs) diff --git a/melange/nn.py b/melange/nn.py index 55db471..cae3ba0 100644 --- a/melange/nn.py +++ b/melange/nn.py @@ -75,7 +75,9 @@ def simple_nnpotential_generator(base_potential, nn_potential, base_potential_pa shapes = [(w.shape, b.shape) for w,b in nn_potential_param_list] def out_nnpotential(x, parameter): - base_e_param, nn_param = parameter[:base_potential_parameter_length], parameter[base_potential_parameter_length:] + base_scale_param, nn_scale_param=parameter[:2] + aux_params = parameter[2:] + base_e_param, nn_param = aux_params[:base_potential_parameter_length], aux_params[base_potential_parameter_length:] #resolve the base param base_e = base_potential(x, base_e_param) @@ -93,11 +95,11 @@ def out_nnpotential(x, parameter): #call the nn_potential nn_e = nn_potential(x, nn_param_list) - return nn_e + base_e + return nn_scale_param*nn_e + base_scale_param*base_e return out_nnpotential -def sigmoid_nn_potential(pos, parameter): +def convolution_nn_potential(pos, parameter, activation=jax.nn.sigmoid): """ make a simple, fully-connected arguments @@ -118,5 +120,5 @@ def sigmoid_nn_potential(pos, parameter): w,b = parameter[i] dotter = jnp.dot(w, act) outs = dotter + b - act = jax.nn.sigmoid(outs) - 0.5 + act = activation(outs) return act.sum() diff --git a/melange/notebooks/csmc.ipynb b/melange/notebooks/csmc.ipynb new file mode 100644 index 0000000..50b3988 --- /dev/null +++ b/melange/notebooks/csmc.ipynb @@ -0,0 +1,825 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "cSMC" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import jax.numpy as jnp\n", + "from jax import random\n", + "from jax import grad, vmap, jit\n", + "from matplotlib import pyplot as plt\n", + "import seaborn as sns\n", + "from jax.config import config; config.update(\"jax_enable_x64\", True)\n", + "from jax.scipy.stats import multivariate_normal\n", + "from melange.gaussians import unnormalized_Normal_logp\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "as a test, i am going to anneal between two gaussians centered at zero over 10 steps. the prior has a variance of 2 and the posterior has a variance of 1.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "DEBUG:absl:Initializing backend 'interpreter'\n", + "DEBUG:absl:Initializing backend 'cpu'\n", + "DEBUG:absl:Initializing backend 'gpu'\n", + "DEBUG:absl:Error initializing backend 'gpu': Not found: Could not find registered platform with name: \"cuda\". Available platform names are: Interpreter Host\n", + "/home/dominic/anaconda3/envs/openmm/lib/python3.6/site-packages/jax/lib/xla_bridge.py:130: UserWarning: No GPU/TPU found, falling back to CPU.\n", + " warnings.warn('No GPU/TPU found, falling back to CPU.')\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "DEBUG:absl:Compiling _where for args (ShapedArray(bool[1,1]), ShapedArray(float64[1]), ShapedArray(float64[1])).\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "DEBUG:csmc:instantiating super\n", + "DEBUG:csmc:define A, b functions\n", + "DEBUG:csmc:getting the parameter lengths\n", + "DEBUG:csmc:getting x dimension\n", + "DEBUG:csmc:getting the terminal time\n", + "DEBUG:csmc:getting the maximum number of twists\n", + "DEBUG:csmc:getting prop_params\n", + "DEBUG:csmc:getting init_params\n", + "DEBUG:csmc:getting model_params\n", + "DEBUG:csmc:getting y\n", + "DEBUG:csmc:define twisting parameters and parameter caches\n", + "DEBUG:csmc:define psi twisting functions\n" + ] + } + ], + "source": [ + "from melange.tests.utils import get_nondefault_potential_initializer, checker_function\n", + "from melange.smc import vSIS_lower_bound, SIS, generate_trajs, SIS_logW, bootstrap_logZ\n", + "import tqdm\n", + "from jax.lax import cond\n", + "import numpy as np\n", + "from melange.csmc import *\n", + "\n", + "#define some parameters\n", + "N,T,Dx = 1000, 10, 1\n", + "potential, (mu, cov), dG = get_nondefault_potential_initializer(Dx)\n", + "cov = cov[0]\n", + "dt=1e-2\n", + "\n", + "#define some _vanilla_ A, b fns\n", + "A_fn = lambda x, params: params\n", + "b_fn = lambda x, params: jnp.zeros(Dx)\n", + "\n", + "potential_params = jnp.linspace(0,1,T)[..., jnp.newaxis] #define the potential parameters\n", + "\n", + "prop_params = {\n", + " 'potential_params': potential_params,\n", + " 'forward_potential_params': potential_params,\n", + " 'backward_potential_params': potential_params[1:],\n", + " 'dt': dt\n", + " } #potential parameters dictionary goes here\n", + "init_params = {'mus': jnp.array([mu]), \n", + " 'covs': jnp.array([cov]),\n", + " 'mixture_weights': jnp.array([1.])} #define the initial params\n", + "model_params=None #there are no model parameters\n", + "y=jnp.arange(T) #the y data aren't important. just know that there are the same number of them as there are data\n", + "rs = random.PRNGKey(6) #define a random number \n", + "\n", + "#create the csmc container object\n", + "smc_obj= StaticULAControlledSMC(N, \n", + " potential=potential,\n", + " forward_potential=potential, \n", + " backward_potential=potential, \n", + " A_fn = A_fn, #A_fn=lambda x, p: jnp.zeros((Dx, Dx)), \n", + " b_fn = b_fn, #b_fn = lambda x, p: jnp.zeros(Dx),\n", + " params_len=1,\n", + " Dx=1,\n", + " T=T,\n", + " max_twists=10,\n", + " prop_params=prop_params,\n", + " init_params=init_params,\n", + " model_params=model_params,\n", + " y=None\n", + " ) #create a controlled smc class\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "can we wrap all of this functionality into something easy?\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "class CSMCWrapper(object):\n", + " \"\"\"\n", + " container for monitoring csmc\n", + " \"\"\"\n", + " def __init__(self, csmc_object, rs):\n", + " \"\"\"\n", + " wrap\n", + " \"\"\"\n", + " assert csmc_object.twisting_iteration == 1\n", + " self.csmc_object = csmc_object\n", + " self.rs = rs\n", + " \n", + " self.Xs_container = []\n", + " self.agg_works = []\n", + " self.bootstrap_logZs = []\n", + " self.out_containers = []\n", + " \n", + " def run_twist_iteration(self):\n", + " re_rs, run_rs = random.split(self.rs)\n", + " self.rs = re_rs\n", + " Xs, works = twisted_smc(self.csmc_object, run_rs, aggregate_works=False)\n", + " bootstrap_logZs = bootstrap_logZ(works)\n", + " adp_fns_dict = self.csmc_object.get_ADP_fn()\n", + " out_container = adp_fns_dict['ADP'](Xs, works, verbose=True)\n", + " self.csmc_object.update_cSMC(new_A_params=out_container['out_A_params'],\n", + " new_b_params=out_container['out_b_params'],\n", + " new_A0=out_container['out_A0'],\n", + " new_b0=out_container['out_b0']\n", + " )\n", + " \n", + " #report\n", + " self.Xs_container.append(np.asarray(Xs))\n", + " self.agg_works.append(np.asarray(jnp.cumsum(works, axis=0)))\n", + " self.bootstrap_logZs.append(np.asarray(bootstrap_logZs))\n", + " self.out_containers.append(out_container)\n", + " \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "DEBUG:absl:Compiling _split for args (ShapedArray(uint32[2]),).\n", + "DEBUG:csmc:generating SMC variable initializer function\n", + "DEBUG:absl:Compiling _split for args (ShapedArray(uint32[2]),).\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "DEBUG:absl:Compiling _where for args (ShapedArray(bool[]), ShapedArray(float64[], weak_type=True), ShapedArray(float64[], weak_type=True)).\n", + "DEBUG:absl:Compiling _where for args (ShapedArray(bool[]), ShapedArray(float64[], weak_type=True), ShapedArray(float64[])).\n", + "DEBUG:absl:Compiling _split for args (ShapedArray(uint32[1000,2]),).\n", + "DEBUG:absl:Compiling _cumulative_reduction for args (ShapedArray(float64[1]),).\n", + "DEBUG:absl:Compiling _uniform for args (ShapedArray(uint32[1000,2]), ShapedArray(float64[], weak_type=True), ShapedArray(float64[], weak_type=True)).\n", + "DEBUG:absl:Compiling _searchsorted for args (ShapedArray(float64[1]), ShapedArray(float64[1000])).\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "Level 9:absl:get_compile_options: num_replicas=1 num_partitions=1 device_assignment=(0,)\n", + "DEBUG:absl:Compiling _normal for args (ShapedArray(uint32[1000,2]),).\n", + "DEBUG:csmc:generating simulation propagation function\n", + " 0%| | 0/9 [00:00" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "sns.distplot(controller.bootstrap_logZs[0], label=f\"uncontrolled SMC\")\n", + "sns.distplot(controller.bootstrap_logZs[1], label=f\"twist 1\")\n", + "sns.distplot(controller.bootstrap_logZs[2], label=f\"twist 2\")\n", + "sns.distplot(controller.bootstrap_logZs[3], label = f\"twist 3\")\n", + "plt.vlines(dG, 0,1000, label=f\"true logZ\")\n", + "plt.ylim(0, 1000)\n", + "plt.legend()\n", + "plt.ylabel('counts')\n", + "plt.xlabel(f\"logZ\")\n", + "plt.title(f\"bootstrapped logZs\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "can we plot the starting and final distributions?" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.stats import multivariate_normal" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "support = np.linspace(-4,4,100)[..., np.newaxis]\n", + "init_untwist_dist = np.exp(multivariate_normal.logpdf(support, mean = mu, cov = cov))\n", + "_, _, (twist_mu0, twist_cov0) = get_twisted_gmm(controller.csmc_object.init_params['mixture_weights'], \n", + " controller.csmc_object.init_params['mus'],\n", + " controller.csmc_object.init_params['covs'],\n", + " controller.csmc_object.A0,\n", + " controller.csmc_object.b0)\n", + "init_twist_dist = np.exp(multivariate_normal.logpdf(support, mean = twist_mu0[0], cov = twist_cov0[0]))\n", + "final_dist = np.exp(multivariate_normal.logpdf(support, mean = mu, cov = np.array([1.])))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXxU1fn/389smewBQsKSAAl72AKyyA5WUCtuVetCLdb6tXbx26/+6tKvrdr22xZbW61Va62t2MVqpdXiVkUFWVVA2dcAgWyQjezLbOf3x72ZTJJJMgkZQpLzfr3yyj1zzzn3uTPJPPc855zPI0opNBqNRqNpjqW7DdBoNBrN+Yl2EBqNRqMJinYQGo1GowmKdhAajUajCYp2EBqNRqMJinYQGo1GowmKdhC9BBHZJyKLuqKuiLwjIitC7CtbRC4OzcrOIyKLRCS3jfNzReSIiFSJyNUduYezuW5PJPAzE5H/FZHnu7DvKhFJN49Xicj/dWHfz4rID7uqP0372LrbAE3XoJSa0Jm6IvIIMEop9ZWA85d1rXXnhB8DTymlfmOWXz8XFxURBYxWSmWdi+sFuf564K9KqU59ySulftaV11FKxXTGjiDXuxW4XSk1L6DvO7uib03o6BGEprcwHNjX3UZ0BDHoFf+DIqIfNnshveKPU9MibPCIiPxDRP4sIpVmSGl687oicinwv8ANZmhgl3l+vYjcbh6PFJEPRaRERIpF5G8ikhCiTZeLyOciUiEiOeZopeHcCBFRIrJCRE6afT8YcD7SDFGcEZH9wIw2rnMUSAfeMO8jotk93Coim0TkMbO/4yJyWUD7r4nIAfO9OiYi3wjx/jaYh7vM694gIv1E5E0RKTKv9aaIpAS0WS8iPxWRzUANkC4iaSKywbz++yLytIj8NaDNhSKyRUTKRGRXQ3hQRH4KzAeeMq//VCt23iIiJ8zP8MFm5x5puJaIOEXkr2a9MhHZJiLJrV3H/Py+LSJHgCMBr40KuESiiKw17+0jERlu1mv4/G0BtqwXkdtFZDzwLDDbvF6Zeb5JyEpE/ktEskSkVETWiMiQgHNKRO4UI+x4xnxPxTw3yrSl3Py7e6X9T7uPopTSP73gB8gGLjaPHwHqgC8CVuDnwMdt1P1rs77WYwzvAUYBS4AIYCCwAXgiWF9BbFoETMJ4EJkMnAauNs+NABTwByASmALUA+PN8yuBjUB/IBXYC+SGcv9B7uFWwA38l/l+fBPIB8Q8fzkwEhBgIcYX97SAe2jrugojRNdQHgBcC0QBscCrwOvN7DoJTMAI8dqBrcBjgAOYB1Q0fCbAUKDE/Cwt5mdRAgxsfp+t2JcBVAELzM/w14An2OcPfAN4w7TdClwAxLV2HfPe15qfUWTz9wNYBVQGXPs3wKZmn7+tjc9sU7PrrQL+zzy+CCgGppl9/xbY0My2N4EEYBhQBFxqnvs78KD5fjqBed39/3u+/ugRRO9lk1LqbaWUF/gLxhdwh1FKZSml1iql6pVSRRhfMAtDbLteKbVHKeVTSu3G+Mds3vZHSqlapdQuYFeAnV8GfqqUKlVK5QBPdsb+AE4opf5gvh8vAoOBZNPOt5RSR5XBR8B7GE/MHUYpVaKU+qdSqkYpVQn8lJb3vEoptU8p5THtmAE8pJRyKaU2AWsC6n4FeNv8LH1KqbXAdgyHEQrXAW8qpTYopeqBHwK+Vuq6MRzcKKWUVym1QylV0U7/Pzc/o9pWzr8VcO0HMUYFqSHa3hbLgT8ppT4z+/6+2feIgDorlVJlSqmTwDog03zdjRGSHKKUqjPfc00QtIPovZwKOK4BnNKJOLGIJInIyyKSJyIVwF+BxBDbzhKRdWa4pRy4M0jb5nY2THIOAXICzp3oqO2tXUcpVWMexph2XiYiH5uhijKML9+Q7rE5IhIlIr83QzoVGCOuBBGxBlQLvK8hQGmATc3PDweuN0M+ZaZ98zAcSyg0eR+VUtUYI5Bg/AV4F3hZRPJF5BciYm+n/5xQzyulqoBS06azZQgBfxNm3yUYI64GWvvbug9jtPipGOHX27rAnl6JdhCa9uR8f27WmayUisN4opUQ+34J42k4VSkVjxFXDrVtAUZoqYFhIbbrECISAfwTI8STrJRKAN4mdDub8/+AscAs8/1a0HCpgDqB73kB0F9EogJeC7zvHOAvSqmEgJ9opdTKIH0Fo8n7aF5nQLCKSim3UupHSqkMYA6wDPhqO9dp7/qB147BCEflA9Xmy4H3PagD/eZjOM+GvqMx7iuvnXYopU4ppf5LKTUEI6z2TLN5E42JdhCa08AIaX01TSxGDLtMRIYC93ag71iMp+M6EZkJ3NyBtv8Avm9O+qYAd3WgbUdwYMSwiwCPOXm9tAPtT2NMkDcQC9RivF/9gYfbaqyUOoERMnpERBwiMhu4IqDKX4ErROQSEbGaE8mLAia+m1+/OauBZSIyT0QcGMuBg37WIrJYRCaZo50KjFCMN8TrtMYXA679E+ATpVSOGa7MA75i3tdtGPNADZwGUsx2wXgJ+JqIZJpO/mdm39ntGSQi1we8f2cwnJG3jSZ9Fu0gNK+av0tE5LMg53+EMRFYDrwF/KsDfX8L+LGIVAIPYXzph8qPMEIIxzHmBP7SgbYhY84T/DeGbWcwnNiaNhs15RHgRTP882XgCYxJ92LgY+A/IfSxHJiNESL5P+AVjAl7zPmXqzBWmxVhjCjupfF/9zfAdeZKnRbzNEqpfcC3Mb5QC8x7bG3j3yAMh1IBHAA+wnBQ7V6nDV7CcJKlGJPeywPO/Zd5LyUYk/ZbAs59iLFs+ZSIFAe5rw8w5lP+ad7XSODGEG2aAXwiIlUYn/V3lVLHO3BPfYaGVRwajeY8wVx2eVAp1eboQ6MJN3oEodF0MyIyQ4z9JhYx9qZcxTnaCa7RtIXe/ajRdD+DMEJ3AzDCP99USn3evSZpNDrEpNFoNJpW0CEmjUaj0QSlV4WYEhMT1YgRI7rbDI1Go+kx7Nixo1gpNTDYuV7lIEaMGMH27du72wyNRqPpMYhIqyoFOsSk0Wg0mqBoB6HRaDSaoGgHodFoNJqgaAeh0Wg0mqCE1UGIyKUicsjM+vRAG/VmiIhXRK7raFuNRqPRhIewOQhTEfJp4DKMrFY3iUhGK/UexdCh71BbjUaj0YSPcC5znQlkKaWOAYjIyxgaM/ub1bsLQ5FxRifaajTnLW6vm6LaInzKSOBmFSuJUYnYLe3l4NFozg/C6SCG0jTbVC4wK7CCmV/gGoz8soEOot22Gs35SmldKXuL93L4zGFcXleTcw6rg/H9xzMxcSLxEfHdZKFGExrhdBDBMnI1F356ArhfKeUVaVI9lLZGRZE7gDsAhg0LS9IxjSYkXF4X63LWcbTsqPlCDdSV4//TFcHlTGBX0S52Fe1ibP+xLEhZoEcUmvOWcDqIXJqmTkzBSBMYyHSM/Ldg5AD+ooh4QmwLgFLqOeA5gOnTp2vlQU23UF5fzjvH36G0phhqiqHyFNSVE2+NIMZM61zpc1Hhy4LIBIgdzCHlo7SulMtGXEaMI6adK2g0555wOohtwGgRScNILXgjzVJOKqXSGo5FZBXwplLqdRGxtddWozlfyK/K5z/Z/6GuphQKD4C7hjR7HBNj0kmxRtMwOlZKkeOpYrermJOFB8ARTVHSeP555J9cOuJSkqOTu/lONJqmhM1BKKU8IvIdjNVJVuBPSql9InKnef7ZjrYNl60aTWc5XX2aN46+gbe6EIoOYVE+FkelMtbRD8QCA0ZCRBwAUlfOsNKjDLPHsre+hI11+aiCnVQPHMeao2u4ZvQ1JEYmdvMdaTSN9Kp8ENOnT1darE9zrqh2V7P68Gqqiw5C2QkixcplUcMZFBEPw2bD4CkQEdu0UV055O+EnE/IdZXzbs1J6pUX+qcTmziW68ZcR6QtsntuSNMnEZEdSqnpwc7pndQaTSfw+Dy8c/wdqkuOwJlsIsTKNdEjGRQzBKZ+FUbMa+kcAJzxkL4QMpeTEj2Iq6LTsGOBkqNUlmbxbva7eH3ec39DGk0QtIPQaDqIUoqPcj6isOQQlB5DRFgaNYyEAaPhghUQG8JcQvxQuOBWEvuN5AtRKcZrxVnkFx9gc/7m8N6ARhMi2kFoNB0kqyyLQ6c/h8KDoBRznYNJ7TcKJn8ZHNGhdxQRC1NuJD1hJDOdg0D5oPAAewu2k12eHTb7NZpQ0Q5Co+kAdZ46Np1cD4X7wedhrKMfk2JHwKTrwNqJ/Qy2CJh0HRfEDCPdHg9eNxTuY0POOtxed1ebr9F0CO0gNJoO8HHBx9QWHQB3LdEWO/NjhiGTrw8+3xAqkf2QSdexMHoYTosNXDVUFe5n26ltXWe4RtMJtIPQaEKkoKqA/blboLIAgPmRQ3CMuwJiB5195wmpRI69jDnOwUa5Io9deZspqik6+741mk6iHYRGEwJen5f1Jz+AkiwARtjjSBt0ASR1ocjw4EzGJk1hqC0GlEIVHWL9yQ/9Yn8azblGOwiNJgT2lezjzKld4K7FLlbmx45AxlwCEkw2rJOIIGMvY2HMMCwi4Kqm6NRODp853HXX0Gg6gHYQGk07uH1uPsvZABV5AMxwJhE76hJwxnX9xaL6kzByKdMikoxy2Um252zQeyM03YJ2EBpNO+wt3ktN0QFQimiLnYlJmTBkavgumDKDzAETjAlr5aOicC8HSw+G73oaTStoB6HRtIHL6+LzE+ugphSA6c5kbGMu7drQUnMsFhxjL2NqxECjXF3M9pPrcPv0slfNuUU7CI2mDXYV7qKu+BAAcRYH41LmQkxS+C8cn8LEIbOIsthBKaqLDrC/RCdU1JxbtIPQaFqhzlPHrpMfQl0FANMjB2FNX3jOrm9Pv4hpTtMZ1ZTyWfaHevOc5pyiHYRG0wp7ivbgKj4CQII1gjHDF0Nkv3NnQMxAMlLnEWNxAFBbfJB9JVr1XnPu0A5CowmCx+dh74kPwVUNwAVRQ7CMmHfO7bClLeKCSHMjXl0Fe06s0/siNOcM7SA0miAcKT1MbamxKS7aYmdU2sUQ0Q1pQSMTGDN8sbGiCagsPsSxsmPn3g5NnySsDkJELhWRQyKSJSIPBDl/lYjsFpGdIrJdROYFnMsWkT0N58Jpp0YTiFKKXSfXQX0VAFOcyVhTZ3WbPfbhs5kYYWaaqytnV85H3WaLpm8RNgchIlbgaeAyIAO4SUSa6xJ8AExRSmUCtwHPNzu/WCmV2Vq2I40mHORU5lBaaMT67WJl/PCFHZPx7mqc8UwcOsfYXQ2cPrWTU9Wnus8eTZ8hnCOImUCWUuqYUsoFvAxcFVhBKVWlGnOeRgO9J/+ppseyK2cj1J4BYHxEfyKGzelmiyBqxALG2M0J8poSdubppEKa8BNOBzEUyAko55qvNUFErhGRg8BbGKOIBhTwnojsEJE7WruIiNxhhqe2FxVp5UvN2VFSW0JO/qcAiAiTBs+CqP7dbBUQm8yU5GnGsVIcz9lCeX1599qk6fWE00EE22raYoSglHpNKTUOuBr4ScCpuUqpaRghqm+LyIJgF1FKPaeUmq6Umj5w4MCusFvTh9lb8AlUGw8aabY44tMXd7NFjQxIv4hUu5F3QlWdYt+pz7rZIk1vJ5wOIhdIDSinAPmtVVZKbQBGikiiWc43fxcCr2GErDSasOH2ujl8ciOYUc/JAydB3JButiqAfiOYlDDGOPb5OHhyPR6fp3tt0vRqwukgtgGjRSRNRBzAjcCawAoiMkrEmHkTkWmAAygRkWgRiTVfjwaWAnvDaKtGw5GSg7hNxdYEawSD0y7uZouaIcKwkUv9G+fqyk6QrZe8asJI2ByEUsoDfAd4FzgA/EMptU9E7hSRO81q1wJ7RWQnxoqnG8xJ62Rgk4jsAj4F3lJK/Sdctmo0APtPrjNyQgMZMalI4uhutqgllqQMxkeZG+c89ew/qZe8asKHLZydK6XeBt5u9tqzAcePAo8GaXcMmBJO2zSaQIpqiigsMpa2WsXC2OEXgeU83EdqtTF+2EK2H/grSilyT++kvL6c+Ij47rZM0ws5D/8DNJpzz/68rX5RvpGOBCJTZnSzRa0Tk3ohw23GZDW1Z9hnrrrSaLoa7SA0fR63182R3MZ9BRmDpnePrEaoRPUnY+Bk41gpDuV8pDPOacKCdhCaPs+Rkv24AienRyzqXoNCYFjaF4i22AGoLTvJ8bKsbrZI0xvRDkLT5zl0Yh2YT+AZcWlIvxHda1AIWBLHkBGdYhS8bg6eWN+t9mh6J9pBaPo05fXlFBQZmdpEhDEjLgpvOtGuwmJhzPDG5EU5hbupcdd0o0Ga3oh2EJo+zeH8bVBfCcAwezxRQy7oZotCJz71QgbbjbkSVXuGw6d3drNFmt6GdhCaPotSikM5G/3lsUmTwRHVjRZ1EGc8Y/ubAslKcVjvidB0MdpBaPosp6ryqSjPBsAhVkYMP390l0Jl5PCFWMX4Ny4uPUJxTXE3W6TpTWgHoemzHM7ZCB4XAKOikrCdhzun2yNi0CTSIky1WXcth/O2dK9Bml6FdhCaPonH5+FIwAazsUNng8XajRZ1EqudsYMbN/Udyd2ic1ZrugztIDR9khOlh3BVnQYgzuJg0LCgavI9gtQRi4gy90RUV+SSZ4bNNJqzRTsITZ/k8ImPwHzSHpOQjsQmd7NFnccSn8roGHNPhM/LIb0nQtNFaAeh6XPUe+s5UbjbXx4zbGEbtXsAIoxJne8vZp/eidvn7kaDNL0F7SA0fY5jp3fiqzPSdQ60RZMw9PwV5guVxJQLibdGAOCqKeZkyaFutkjTG9AOQtPnyMppFOYbPWA8OKK70ZquQaL6MbqfuQpLKY6c3Nh2A40mBLSD0PQpalzV5Bbv95dHBoRmejqjUub6j08U7cbldXWjNZreQFgdhIhcKiKHRCRLRB4Icv4qEdktIjtFZLuIzAu1rUbTGY6d2oFyVQMw2BFL7ODMbrao6+g/dCYDbMZOcG9tGccD5lk0ms4QNgchIlaMNKKXARnATSKS0azaB8AUpVQmcBvwfAfaajQdJjDvw6iBU8Dm6EZruhhHFKMGjPcXs3SYSXOWhHMEMRPIUkodU0q5gJeBqwIrKKWqzBzUANGACrWtRtNRquorOVV6GDCUW0f24L0PrTFqWGPILKd4P3Xu2m60RtPTCaeDGArkBJRzzdeaICLXiMhB4C2MUUTIbc32d5jhqe1FRUVdYrimd5KVtwXlrgNgaEQCUUm9b1AaP2gqyQ4jHanPVcWxgh3dbJGmJxNOBxFMVF+1eEGp15RS44CrgZ90pK3Z/jml1HSl1PSBAwd22lhN7ycr92P/8ajkaT1TWqM9bA5GNaQjBbJyN3WjMZqeTjgdRC6QGlBOAfJbq6yU2gCMFJHEjrbVaNqjoq6MwrKjgBFeSh/ewzfHtcHI1MbQWV7pYWp1IiFNJwmng9gGjBaRNBFxADcCawIriMgoESN9l4hMAxxASShtNZqOcCx3C3iN3cUpzkSc/Ud1s0XhIyZ5AoMc8QAody3HA0QJNZqOEDYHoZTyAN8B3gUOAP9QSu0TkTtF5E6z2rXAXhHZibFq6QZlELRtuGzV9H6OBXxJjkye2jPSinYWi5X0pMYw09G8j9uorNG0ji2cnSul3gbebvbaswHHjwKPhtpWo+kMVXXlnAoIL6UN6z2b41pjZOo8tuQay1xzSw5R567FaY/sZqs0PQ29k1rT6zmauxm8HgBSIgcS2X9kN1sUfmKTJpAcYYaZPHUcz/+kmy3S9ETCOoLQaM4HAsNL6cmZYQ8v+XyKyjoPZ2pcVNV7aNjpIwKxThsJUQ7inDYknHZYrIwcOJnT5ijiaN4njB++KHzX0/RKtIPQ9Gqq6yuahZe6fvWSUorSahfHi6s5VlzN6fI6PL6gq7L92CzC4IRI0hKjSU+Mpl901+/oTh82PyDMdFCHmTQdRjsITa/mWM4mlBleGuJMJKpfWpf17fb6OHSqks9zyiiurO9QW49PkVNaQ05pDRsOF5EUF8HU1H6MSY7BZu2ayG/cwAkMjIinqL4cn6eO7PyPGTd8cZf0rekbaAeh6dUcDQgvjRo0rUvCSy6Pj89OnmFnThm1Lm/QOtERVhKiHMRH2rGa1/T4FBW1bs7UuKhp1q6wop53951iU5aVqcP6MTU14ewdhcXCqKQpFOVsAOBo7ifaQWg6hHYQml5LTX0lBU3CS2envaSU4tDpSjYdKaayztPknN0qjEiMJi0xmhEDoomOaPtfq7LOzYmSGo4VV3OiuNofkqqu97LpSDF7cstZMCaRkQNjzmquIj11HltNB5FTehCXpx6HLaLT/Wn6FtpBaHotx3O3+MNLg84yvFRe4+bd/afIO9NU/C7WaWPqsAQmDInHaQ9duiPWaWfi0HgmDo2n1uVlT145u3LKqKo37C2vdfPGrgKG9Y9iyYRk4pz2TtkdP3ACiRHxFNeX4/PUcyL/U0b3gWW+mq5BL3PV9FqOFWzzH6cnT+lUeEkpxf78Cv76yYkmziHKYWVJRjK3zU3jguH9O+QcmhPpsDIzrT+3zUvjonFJTfo6WVrD3z4+yZHTlZ3r3GIhbeBEf/FYnl7uqgkdPYLQ9Erq3TXklR7xl9M7kTnO5fHx/oHTHDrV+OVsEWHqsARmpfcnwta1Yn9WizAlNYGxg2LZerSEXbllKAV1bi9v7i5gwpBqFo9Lwt7BuYn0lDlsM/NgnCzZj9vrwm7tRXkwNGFDjyA0vZLsvI/xmSk3Bzr7EztgdIfaV9S5+cf2nCbOoV+UnRtnprJgzMAudw6BOO1WFo9L4vrpqcRFNoaW9uVX8M8duVTXe9po3ZL+SZOJd8QA4HbXklvwWZfaq+m9aAeh6ZUcDwilpCdN7lB4qaC8lpc/PUlRwNLViUPjuXnWcJLjnF1qZ1sMTYhk+axhjBsUG2BbHX9vZlt7iNVGeuIEf/lY3tYutVPTe9EOQtPrcLvrOVl6yF8emTqvjdpNOVpUxertuVTXG8tQLSIsyUhmSUYyDtu5/3dx2q1cOnEQi8YO9Pu4yjoP/9iew8mS0GW8R6bM9R9nF+/H6+3YKETTN9EOQtPrOFnwCR6P8YTdPyKehIGhZY47fLqSN3cV+JecOu1WvjRtKBOHxofN1lAQEaYO68fVmUP9Tsrl8fHvnXkcK6oKqY+Bg6YSY4sCoN5VRf7pXWGzV9N70A5C0+sIXKmTPjC08NK+/HLe3lOAzxROio+0c9PMVFL7R4XNzo4yIjGaG2ekEus01pZ4fIo3dxeEtMLJCDM1OspjeVvCZqem96AdhKZX4fV6OFF8wF9OT5nTbpt9+eW8t++0X1Svf7SD66enkBB1/q30GRATwfUXpBJvTl57fYq39oTmJNKHXug/Pl60F+Xzhc1OTe9AOwhNryKvYAcuj7FfIc4Rw4CkSW3WzyqsZO3+0/5yYmwE109PIbaTG9POBfFRdq6fnkK/KMNGpeCdvafILq5us92gITOItBq7qGvqKzhdpHNwadomrA5CRC4VkUMikiUiDwQ5v1xEdps/W0RkSsC5bBHZIyI7RWR7OO3U9B4CV+ikJU5ArK0vRz1ZUsPbe075Rw4DYyO4/oIUohzn//agWKed66en0t9UgfX6FG/uzie/rLbVNhabgxEDxvrLx3J1mEnTNmFzECJixUgjehmQAdwkIs1nC48DC5VSk4GfAM81O79YKZWplJoeLjs1vQefz8vx4san4vShs1ute6q8jjd25+M1J6T7Rdm5ZurQs9oRfa6JjrBxzbSh/jkJt1fx+s68NpfApg9pDDMdK9qNUm3Lkmv6NuEcQcwEspRSx5RSLuBl4KrACkqpLUqpM2bxYyAljPZoejmnC/dQ6zJW9UTZI0kePC1ovfJaN//emYfLY8TgY502vnRBSrsCe+cjcU47105LIcphOLZ6t7G6qaqVzXRDh87CYTFCUxW1pZSUZp0zWzU9j3A6iKFATkA513ytNb4OvBNQVsB7IrJDRO5orZGI3CEi20Vke1FR0VkZrOnZBIZMRvQfj8Xach6hzu3l3zvz/HLbkQ4rX5qW0mkxvPOBftEOrpnWuAS2ss7Dmp35fgcYiM0RxbB+o/zl4zmbzpmdmp6HhGuIKSLXA5copW43y7cAM5VSdwWpuxh4BpinlCoxXxuilMoXkSRgLXCXUmpDW9ecPn262r5dT1f0RVz19Xyy8z84bP0QINIRi83WdNezUlDr9uDxNuQAhSiHDZslvClIzxUer48at9d4tAJsVot/ZBGI211LndsYaVnESnRk/3NppqabcDqdpKSkYLc3fRgSkR2thfHDOabOBVIDyilAfvNKIjIZeB64rME5ACil8s3fhSLyGkbIqk0Hoem7HM8+zJDkUcTEObFYLCTGDUcsTQfIFbVuat2NiXriIu1E9qA5h1CocXma5KqIclhbrMjy+byUVJxoXNYbOxSbzhHRq1FKUVJSQm5uLmlpocvehzPEtA0YLSJpIuIAbgTWBFYQkWHAv4BblFKHA16PFpHYhmNgKbA3jLZqeji1dbXExDkRESKszhbOodblaeIcoh22XuccwBgRBY4aalxe6txNs9dZLFbs1kaH4HKHthtb03MREQYMGEBdXV2H2oVtBKGU8ojId4B3ASvwJ6XUPhG50zz/LPAQMAB4xsya5TGHOsnAa+ZrNuAlpdR/wmWrpufj8br9mdcc9qa7n91eHxUBT9URNgvREb3POTQQE2HD41P+OYiKWjc2izRJYRphi8JlypHUu2uIihzQLbZqzh2dyUwY1mUbSqm3gbebvfZswPHtwO1B2h0DpjR/XaMJRnlZNj5lPCWLgMOUtgZjf0BZjdtftlmEuEj7WaXxPN8REeKddkprXHh9CgWU1brpH+XAYs63RDhiqawzFhC6vW68XjfWIJP6mr6N3kmt6fEcO9k4NeWwRmCxGKMDpRQVtW6/vpKIobFkCYNzKCsr45lnnulQm+PHjzNr1ixGjx7NDTfcgMvl6hJbXn/9dQ4ePEB8pJ2GO/X6FBV1bv++B4vV5k8a9OILL/HCC39otb/169ezZUvHN9WNGDGC4uLiNuusWrWK73znOwA8++yz/PnPf+60HWvWrGHlyiuxVvAAACAASURBVJUA3HrrraxevbpD9v7sZz9rUp4zp32Zlt6OdhCaHs/x0zv9xw57tP+42uXF5W1c6hnvtDcJs3QlnXEQ999/P3fffTdHjhyhX79+/PGPf+wSW15//XX279+P3WppknCo3uOj1tU4HxFhhuJWfO1mrvvyla3211kH0VHuvPNOvvrVr3bKDo/Hw5VXXskDD7QQbAiZ5g7iXNzz+Y52EJoeTXVFPqeqCoyCQIQZXnJ5vE0yr0U7rESEcVL6gQce4OjRo2RmZnLvvfe2W18pxYcffsh1110HwIoVK3j99ddb1At8wgZYtmwZ69evByAmJoYHH3yQKVOmcOGFF3L69Gm2bNnCmjVruPfee8nMzGTXZ9u5ZKHxJLxvz26inXaOHs8GYNLE6dTU1PLLlb/hiSeexufz8OSTT5KRkcHkyZO58cYbyc7O5tlnn+Xxxx8nMzOTjRs3UlRUxLXXXsuMGTOYMWMGmzcb6UxLSkpYunQpU6dO5Rvf+Earu7RfeOEFxowZw8KFC/1tAR555BEee+wxgJDsuPXWW7nnnntYvHgx999/f4v36v3332f+/PmMGTOGN998s83384EHHqC2tpbMzEyWL1/uf38bPqt7772XiRMnMmnSJF555RXAcFiLFi3iuuuuY9y4cSxfvrzX7UwPaQ5CRJYBbyultPyj5rzieM5G/7HD4sBisfHr9w5R6/b6l3FaLEKk3QKcfWjp7iVjgr6+cuVK9u7dy86dxmimsrKS+fOD58F+6aWXSEpKIiEhAZvN+BdMSUkhLy+vQ7ZUV1dz4YUX8tOf/pT77ruPP/zhD/zgBz/gyiuvZNmyZX7n43bVU1tVycdbNzNl6jTWfrAeyxcWkZSUTFxsnL8/l6uKlStXcvz4cSIiIigrKyMhIYE777yTmJgYvve97wFw8803c/fddzNv3jxOnjzJJZdcwoEDB/jRj37EvHnzeOihh3jrrbd47rnmyjlQUFDAww8/zI4dO4iPj2fx4sVMnTo16PvZnh1//OMfOXz4MO+//z5Wq5VVq1Y16SM7O5uPPvqIo0ePsnjxYrKyWt81vnLlSp566in/5xfIv/71L3bu3MmuXbsoLi5mxowZLFiwAIDPP/+cffv2MWTIEObOncvmzZuZNy/0BFXnO6FOUt8I/EZE/gm8oJQ60F4DjeZccPxUY35lhz0apRT1Hp/fOSDgtHWNc+gIsbGxQb9sGgi267+jE+cOh4Nly5YBcMEFF7B27dqg9ebMmcPezz7lk82b+O7/u491769FBObNm0eErXHFV727msmTJ7N8+XKuvvpqrr766qD9vf/+++zfv99frqiooLKykg0bNvCvf/0LgMsvv5x+/fq1aPvJJ5+waNEiBg4cCMANN9zA4cOHW9QLxQ6A66+/Hmsrgoxf/vKXsVgsjB49mvT0dA4ePNhqP22xadMmbrrpJqxWK8nJySxcuJBt27YRFxfHzJkzSUkxFIIyMzPJzs7uew5CKfUVEYkDbgJeEBEFvAD8XSnVvhC9RhMG6qqLya1sUHMRIhwx1Lq9fgE+MJxDd6xYam8EMX78eMrKyvB4PNhsNnJzcxkyZEiLujabDV9A3obAdex2e+NqLKvViscTXH9p/vz5bNmymfy8HC69/AqeeuJXiAhXXrHMH5IDcHnreeONNWzatJk1a9bwk5/8hH37WkqC+3w+tm7dSmRkZItzobzXodR566232LBhQ5t2AERHRwd9Pdh1RKTN97M12gobRUQ07idp6zPoqYS8zFUpVWGOICKB/wGuAe4VkSeVUr8Nl4EaTWucyNno/+e1Wqz4sFJZ5+K/FqQDwXcRh4vY2FgqKyublNsaQQAsXryY1atXc+ONN/Liiy9y1VVXtagzYsQInnnmGXw+H3l5eXz66acdtmXBggX84Ac/YMGCBURH2Eno158P1v6HBx/+MYgdEWMq0uv1cezYQRYvXsy8efN46aWXqKqqIjY2loqKCn9/S5cu5amnnvLPtezcuZPMzEwWLFjA3/72N37wgx/wzjvvcObMGZoza9Ysvvvd71JSUkJcXByvvvoqU6Y0XdHu8/nIyclp1472ePXVV1mxYgXHjx/n2LFjjB07lsrKylbfT7vdjtvtbiFFsWDBAn7/+9+zYsUKSktL2bBhA7/85S87PSLpSYQ0SS0iV5pyFx8CdgxNpcsw9ip8L4z2aTStcqxgh//YanFQUdt0v0PMOVRnHTBgAHPnzmXixIkhTVIDPProo/z6179m1KhRlJSU8PWvf71Fnblz55KWlsakSZP43ve+x7RpwRVqA7nxxhv55S9/ydSpUzl69CgjRowAjC+6GKeNC2fPIS4+gfh+/aioc2OzGO+T1+vl1ltvZ9KkSUydOpW7776bhIQErrjiCl577TX/5PCTTz7J9u3bmTx5MhkZGTz7rLG16eGHH2bDhg1MmzaN9957j2HDhrWwbfDgwTzyyCPMnj2biy++OOj9eL1evvKVr7RrR3uMHTuWhQsXctlll/Hss8/idDrbfD/vuOMOf2grkGuuuYbJkyczZcoULrroIn7xi18waNCgdq/fGwhJrE9E/gw8H0wsT0S+oJT6IBzGdRQt1td3cNeW88J738bjM5ZtXpj+AEPTGpPh9I92YA/Tktaejtvjo7Smcc+F0+alvt7IqiciQXWsNL2DAwcOMH78+CavtSXWF+pfQUFz5yAijwKcL85B07c4mbPJ7xwi7Im4vY3x5linTTuHNrDbLE1GV3UeK0Z+LyPe7tbaTBqTUP+LlgR57bKuNESj6QjHCrYB4FOKSnejOqXDZumVInxdTZTD2sSJ+pTDf1zvbju3tabv0KaDEJFvisgeYFxA7ujdInIc2H1uTNRomuKpr+RE2VHAyA7ntmcChpRGnLN36yx1FSJCnNPmX/yrJNK/+qveU4fe8qSB9lcxvYSR5e3nQOAe9kqlVGnYrNJo2iA3ZwsunyHfXeeNI9IxFKgn1mnH2kuS/5wLbFYLMU6bkT9CIvApwaIUPsDtqsYREdvdJmq6mfZCTEoplQ18G6gM+EFEdBoqTbdwLH8bPqU4U+MmMmIsIoLdajE3xGk6QqTdisPcK2KxOPE0jCJ0mElDaCOIZcAOjESGgY9nCkgPk10aTVB8rhqyy45QVuvGqxS26GlER1hx2rtnQ1xPpyHUVFLtQixRKF8NXp+i3lNHjPL590ho+iZtfvpKqWXm7zSlVLr5u+FHOwfNOSc/92PO1NVS7fJiscZjcwznonHJ3e4cwiX3vXPnTt5+++0grYOzfft2/vu//7vV89nZ2bz00ktNXrNaLMRG2EEiAAtepfD6vHjcNf46ochnZ2dnM3HixE7bEUh+fr5fS6q5wF4orFq1ivz8xgzHt99+exN5EE1ohLpRbq6Z+hMR+YqI/NpMF9peu0tF5JCIZIlICx1eEVkeMPG9RUSmhNpW0zc5lLuVMnNDXJRjDOMGxzEqKaadVuEnXHLfHXUQ06dP58knn2z1fGtfzE67hQi7FYsY0hEen6LO1fnlrp21Awzp7iFDhnQ4n0MgzR3E888/T0ZGRqf766uEOn78HVBjfoHfB5wA/tJWAzEWVj+NsRw2A7hJRJp/QseBhUqpycBPgOc60FbTx1Duenbm7fPHyWMSLmDx2KRutsogHHLfLpeLhx56iFdeeYXMzExeeeUVJk2aRFlZGUopBgwY4E+wc8stt/D++++zfv16v4DfRx99RGZmJpmZmUydOpXKykoeeOABNm7cSGZmJo8//jher5d7772XmTNnsvDC6fz1z/8AjOXD3/2f+8jIyODyyy+nsLAw6H3s2LGDKVOmMHv2bJ5++mn/6x21Y9WqVVx//fVcccUVLF26tMloBCAnJ4dLL72UsWPH8qMf/QigRZ3HHnuMRx55hNWrV7N9+3aWL19OZmYmtbW1LFq0iIZNtH//+9+ZNGkSEydO5P777/e3Dyaf3tcJVYvAo5RSInIV8Bul1B9FZEU7bWYCWWb6UETkZeAqwD/OU0oFZuT4GEgJta2m77Hr4EcU1xlhD7HEcGnmbCIdQfY8rPt5+IxY/P2gL4dD7tvhcPDjH/+Y7du389RTTwGwbt06Nm/ezPDhw0lPT2fjxo189atf5eOPP+Z3v/sdgUoCjz32GE8//TRz586lqqoKp9PJypUreeyxx/z5EZ577jni4+PZtm0b9fX1zJ4zh3kLprJvzz6yso7y6batVJRXk5GRwW233dbiXr72ta/x29/+loULF7bqGEOxY9WqVWzdupXdu3fTv39/srOzm/Tx6aefsnfvXqKiopgxYwaXX345iYmJQa933XXX8dRTT/HYY48xfXrTDcL5+fncf//97Nixg379+rF06VJef/11rr766lbl0/syoY4gKkXk+8BXgLfMJ/z2VNCGAjkB5Vzztdb4OsaS2g61FZE7RGS7iGwPJqGs6R24PD4271vnL6cmTmTsoPhutKhtGsT6gv1kZGQEVQgNZR5l/vz5bNiwgQ0bNvDNb36TPXv2kJeXR//+/f0JbhqYO3cu99xzD08++SRlZWV+ZxTIe++9x5///GcyMzOZNWsWZ0pLyTlRwCdbt3HlNcuora9m8ODBXHTRRS3alpeXU1ZWxsKFCwFjFBOMUOwAWLJkCf37B18cuWTJEgYMGEBkZCRf+tKX2LRpU5vvU2ts27bNLzdus9lYvnw5GzYYIhHN5dObO6m+SKgO4gagHvi6UuoUxpf1L9tpE+yvPajwk4gsxnAQDeO9kNsqpZ5TSk1XSk1v0JjX9D62HMqluNZI+GIRuGjKJd1sUdtUVlb6wyrNf/bv309iYqJf7htoVe67OQsWLGDjxo1s3LjR/0W3evXqoKOVBx54gOeff57a2louvPDCoOqjSil++9vf+p3X8ePHueKLxpekiOD21VFd727RrqFtKE4tFDug+6W7Q5VP70uEmg/iFPDrgPJJoPXs4ga5QGpAOQXIb15JRCYDzwOXKaVKOtJW0zcoKK9l3+FN1GF8UQ2MT2TUkMmtN2glDBROwiX33bzf1NRUiouLcblcpKenM2/ePB577DF/CCqQo0ePMmnSJCZNmsTWrVs5ePAgqampTfq75JJL+N3vfsdFF12E3W7n8OHDDB4ymNlzZvHnVS9x3Q3XcDz7KOvWrePmm29u0n9CQgLx8fFs2rSJefPm8be//S3ofYZiR3usXbuW0tJSIiMjef311/nTn/5EcnIyhYWFlJSUEBMTw5tvvsmll14a9H1roEFuvLi4mH79+vH3v/+du+66K2Q7+hqhrmL6kogcEZFyEakQkUoRaU+YfRswWkTSRMSBkZVuTbN+hwH/Am5RSh3uSFtN38Dj9bF2/2lc1bsAIwHQ5NRpWC3nl95SuOS+Fy9ezP79+/2T1GB8yY0ZY6Q+nT9/Pnl5eUGzmD3xxBNMnDiRKVOmEBkZyWWXXcbkyZOx2WxMmTKFxx9/nNtvv52MjAymTZvGxIkT+cY3voHP6+PaL11NevoIliy4nPvuuZvZc+cHffp+4YUX+Pa3v83s2bODJhAK1Y72mDdvHrfccguZmZlce+21TJ8+HbvdzkMPPcSsWbNYtmwZ48aN89e/9dZbufPOO/2T1A0MHjyYn//85yxevJgpU6Ywbdq0oI5ZYxCq3HcWcEVHU42KyBeBJwAr8Cel1E9F5E4ApdSzIvI8cC3GqigwJsOnt9a2vetpue/ex9ajJXx6JJ/ygh9TL26SY51cPf9/GZbUdAQRTMZY03nqXFWUV5/G7VVYsGC1DybGaSf6HObY0HQ9HZX7DvXTPt2ZPNRKqbeBt5u99mzA8e3A7aG21fQtSqrq2ZZdir1mD7W4SXDaiY5KYGjihO42rdfjsEdhEStW8eJVPmy+OqrrhQibBZuWUu8zhOogtovIK8DrGJPVACil/hUWqzR9HqUUHxwoxOtT+Gp2EWE1chikJU0+78JLvRGLWHDYnChVg8+rUL4asEZSWe8hIVIr5vYVQnUQcUANsDTgNYUxf6DRdDm7c8vJK6vF4q2jwn2cfrF2BBiZ0jLWrgkPEY4Y6t012CyCz1ePRSlcHh91bl/w/SeaXkeoq5i+Fm5DNJoGKuvcbMoqBsBRvRub04fdasfhjCdloA4vnSsc9ihTrM8HohBfHcoaSWW9mwibBYuWVu/1hLqKaYyIfCAie83yZBHp21sMNWFj3aEiXB5jfbvFvZs4p7EnMy1pig4vnUOMMJOxMslqEVDGLnaloLJe7xHoC4Q62/QH4PtgLERXSu3GWHqq0XQpWYWVHC00ROKsnhqwF/h3TY4ZvrD7DOujRDgad2eLuMDMNFfn9lLv8XaXWZpzRKgOIkop9Wmz1/QjhKZLqXN7WXewUS5lVMRBvBbjz8wZNYChA8a11rTb6Yya64cffujff7BixYpWd+5+/vnn3H67sdivM9LX0FI9dc+ePdx6662t1r/pppuYPHkyzzz1ex79+RN8tH4zCkWENEqSV9Z68LWzTD5QtC+QjqrUng3NlV1bQ8uVtyRUB1EsIiMx5S5E5DqgIGxWafokm7OKqTJDF9ERVpwB2owjk6diOY+T13TUQfh8PlasWMHLL7/M3r17GT58OC+++GLQuj/72c/Oerdv8y+0SZMmkZuby8mTJ1vUPXXqFFu2bGH37t3cc889PPzwD1i4aC4AVqmlYQGTVymqOxlq6oyD6Kz0RagOIhAtV24Q6n/ct4HfA+NEJA/4H+DOsFml6XPkldWyO7fcX1403MGJyuNGQYRRwxZ0k2Wh0VG575KSEiIiIvw7opcsWcI///nPFvUqKyvZvXs3U6ZMaXHujTfeYNasWUydOpWLL77YL08dirw2wBVXXMHLL7/cot+lS5dSWFhIZmYmGzdu5NvfvIc3/m3oaE6YMIMnVv4fS+bPZtHs6ezaux+318enn37KnDlzmDp1KnPmzOHQoUOt3nswGfPW2jeXAff5fHzrW99iwoQJLFu2jC9+8Yv+L+IdO3awcOFCLrjgAi655BIKCgqCSn8HouXK26bNVUwick9A8W1gHYZTqcbYAf3rYO00mo7g8fr44EDjH3P6wGictdup9ZmjiehkBvcbGXJ/z+zsWKinI3wr81tBX++o3Pf48eNxu91s376d6dOns3r1anJyclrU3b59e5MvkUDmzZvHxx9/jIjw/PPP84tf/IJf/epXIclrg/GUvHLlSu67774m/a5Zs4Zly5b57+WPf/wjYml8lhw4IJ71Wz7huWd/x++efILxzzzL2LFj2bBhAzabjffff5///d//DerwILiMeUVFRavtA2XAV69eTXZ2Nnv27KGwsJDx48dz22234Xa7ueuuu/j3v//NwIEDeeWVV3jwwQf505/+1Kr0N2i58vZob5lrrPl7LDAD+DeG0uotwIazurJGY7It+wwlVUZs22GzsHhcEtu3fOI/P3LQtPM6vBSMUMT6Xn75Ze6++27q6+tZunRpUBnsgoICWlMpzs3N5YYbbqCgoACXy0VaWhrQKK+9fPlyvvSlL5GSkhK0fVJSUsihF5s1wn/8xcsvIs5pY8rUqbz9xr/x+BTFRaU8cOutHDlyxFCBdQdXgG2N8vJyVqxYEbR9oAz4pk2buP7667FYLAwaNIjFixcDcOjQIfbu3cuSJUsA8Hq9DB48uN1rNpcrf+edd1rUC/X9DEWuHPDLlV999dVt2heMQLlywC9XfvXVV7eQK1+7dm2H+29Omw5CKfUjABF5D5imlKo0y48Ar5711TV9ngY5jQbmjBxAlLuUYxVmbFwsjOqBq5faG0FkZGQwe/ZsNm7cCBh5GQ4fPtyibmRkZKsy1nfddRf33HMPV155JevXr+eRRx4BjHDX5Zdfzttvv82FF17I+++/H7R9XV1dqwJ7zXFYnY0FG1jFR7TTgcdrjPIefvghFixaxGuvvUZ2djaLFi0Kqd8GfvjDH7J48eKg7QNlwFvTjlNKMWHCBLZu3RryNTsiVx7K+9kb5cpD3Uk9DAjMqO4CRpz11TV9mkA5DYBB8U6mpCRwYv+ruJSxhDI2djDJsaltddOC1sJA4aQzct+FhYUkJSVRX1/Po48+yoMPPtiizvjx4/nVr34VtH15eTlDhxp5tAInuEOV1z58+HCr4avmWK12LBbj68KnFO76KiLtVv8S5IqKcvoPHIRSilWrVrXbX/P3K/Be2mo/b948XnzxRVasWEFRURHr16/n5ptvZuzYsRQVFbF161Zmz56N2+3m8OHDTJgwoVXpby1X3j6hjtv/AnwqIo+IyMPAJ0DwJRcaTYjsyTPkNAAsIlw8PhmLwOH8xvDSqCGzeoTuT2fkvn/5y18yfvx4Jk+ezBVXXBE0a9u4ceMoLy8P+mXxyCOPcP311zN//vwm8exQ5bXXrVvH5ZdfHpKtIkKEPcpfrnNXIYDNnJv49nfv4ccP/YA5c+fi9ba/P6K5jPl9993H97//fea20/7aa68lJSXFL00+a9Ys4uPjcTgcrF69mvvvv58pU6aQmZnJli1GRuPWpL9By5W3R0hy3wAiMg1oGDNvUEp9HjarOomW++45VNV7eHFLtn/H9My0/swdlUh9SRarNj6MV/nAauOGL/yaAdHtZwrszXLfjz/+OLGxsf69EF1BfX09CxcuZNOmTa2mAG2O2+viTGUuKGMiMjE2BbE5qKh1U+s2vtQtAgOiI8Iqw1FVVUVMTAwlJSXMnDmTzZs3M2jQoLBdrzcRLrlvlFKfAZ+dnXkajcG6g4V+55AQZWdmmjG5dyx7neEcgMSE9JCcQ2/nm9/8Jq++2rVTfidPnmTlypUhOwcAm8WO1eLA63WhgHpXJU7bAGKcNuo9PnxK4TNlOOIj20tZ33mWLVtGWVkZLpeLH/7wh9o5hJGwZv8QkUuB32Ak/XleKbWy2flxwAvANOBBpdRjAeeygUrAS0AiIU3PJ6uwkixTTgPg4vHJ2K0W8Lo5fLrxGWRMqlZuBXA6ndxyyy1d2ufo0aMZPXp0h9qICE5HNNW1xnRkvbsap+qPRYRYp43yWmPlUZ3bi9NuIcIWHt2s9evXh6VfTUvC5iBExAo8DSzByDG9TUTWKKUC94uXAv8NtLbea7FSqjhcNmrOPc3lNCYOjSe1vxHbrjy1m3yXsVlO7JGMHjq7Q32HuipF03kiHLFU15WBUtT7PPg8dVjskTjtVlOfyRj9VdZ6sMdYsOjP47wh1OmEQMK5uHwmkKWUOqaUcgEvA01mU5RShUqpbZgigJrez6YjTeU05o9unFzNOrHe/0eckphBdIBQXHs4nU5KSko69U+gCR2bxYY9YE9EvStg5ZbT3iUyHJquRylFSUkJTqez/coBhDPENBQI3BqaC8zqQHsFvCciCvi9Uuq5rjROc+7JKa1hT16AnMbYJJx2MwzhquZwSePgckwHpTVSUlLIzc2lqKio/cqas8LlqfM7BqtYiHIW0+AZjIRC5iokgWiHzZAK13Q7Tqez1U1+rRFOBxHsr6Ijj3dzlVL5IpIErBWRg0qpFru3ReQO4A6AYcOGdc5STdhxe328HyCnMTIphtFJjSOE4pNbKfEY+QZskQmkJWd2qH+73e7fSawJLzWual5c+12UqxqA5RfcRfywOYDxpPqvz/I4WWp8lomxEdw8c5h2Ej2UcIaYcoHAHU4pQMiSikqpfPN3IfAaRsgqWL3nlFLTlVLTW5Mk0HQ/Hx8roazGiCRG2C1cNC6pcb5AKQ7mrPfXTUueisPq6AYrNaEQ5YhmWNJkf/lQ9of+YzH3s9itxmdbXFnP9oCd8pqeRTgdxDZgtIikiYgDI8HQmlAaiki0iMQ2HGPkwt4bNks1YeV0RR2fnSjzlxeMHkhMROPg1VuRx+GKE0bBYmHciC+caxM1HWTc8MX+sNLBM0dQVY1rSeKj7Mwe2Ti39MnxUkqq6s+5jZqzJ2wOQinlAb4DvAscAP6hlNonIneKyJ0AIjJIRHKBe4AfiEiuiMQBycAmEdkFfAq8pZT6T7hs1YQPr0/x3v7T/sQyqf2jmDAkrkmdE8c/pM5Ubo2JS2Vogg4Vne+MSMzAGZ0MQJXPRe6JD5ucn5qawKB4Y0LU61Os3X8an08vIOhphFUiUyn1tlJqjFJqpFLqp+ZrzyqlnjWPTymlUpRScUqpBPO4wlz5NMX8mdDQVtPz+PR4KcWVxtOj3SpcPD6p6VJUr5sDBY3SGmNT5vY45da+iNViZUzAMuSDuVvB1yiRYbEYoaaGuYeC8jo+zylr0Y/m/Eb/J2rCRlFlPZ8eb4w/zx6ZSEJU07mF6oKdnKw/YxTskYxLDa6Aqjn/GDd8EZhzRcdqC6krOtDk/MDYCGaMaJS/3pJVzJlqF5qeg3YQmrDgawgrmKGlIQlOpqYmtKh3KPtD/96FIYkZxDtb1tGcnyRGJ5HY39iN7VU+so5/0KLOzLT+JMYa+yY8PsXaA6f1XpUehHYQmrCw/cQZTlcY+vY2i7AkY1ALATdVU8rBUvOpU8R4ItX0KMaPWOw/Pli0B+qrmpy3WoRLMpL9O6rzztSyU4eaegzaQWi6nKLKej4+VuIvXzhyAP2jWy5bPX3iI8q85vxEVCIjk1rmXdac34wedAEWc9RX6Kmm+OTmFnWS4pxMH9HPX96cVUxZjQ419QS0g9B0KcaqpVNNkgBdMKxfy4o+L/tONu57HDVkFnZr+BRANeHBaXOSPvgCf3l/9ocQkCWtgVlp/UmMMR4S3F7Fe/v0qqaegHYQmi7l0+OlFFYYowKbRViakRw0N0DdqT1k1RYaBVsEGSNaJsvR9Awy0paA1djXcrg6D3dJy9SpNquFpRMGNYaaymr1qqYegHYQmi6jsKKuyaqlOaMGMCAmImjdg8febcz7MGAsSdFa07+nMjRuGPEJ6QC4lJcjx9YGrZcc52RGWuNocktWMaV6VdN5jXYQmi7B4/Xxn32n/KuWhiZEMjU1SGgJmL3XTQAAIABJREFUUFVF7AuYnJ6QdrGW6e7BiAgT0hp3v+8r3Al15UHrzkobwMCAVU3v7msMR2rOP7SD0HQJm4+WUFJlPA3arcKSVkJLAHnHP6DcnJx2RCcxJqljwnya84+xg2dgjTL2PBR5aig8vj5oPatFuGTCIP8GulPlTUedmvML7SA0Z01OaQ2fnTjjL88fPZB+QVYtAeBxsS+3caXLmJS5enK6FxBpi2TUkEY1/305G5rsrA5kYGwEs0cO8Jc/PV7KqfK6sNuo6TjaQWjOijq3l3f3nfKXRyRGMTklvtX6Nfk7OFZvLoG1R5IxfHGrdTU9iwkjvgA248HgSM0p6k/tabXuBcP6MTQhEgCfUvxnbwFub8vVT5ruRTsIzVmx/lAhlXWG0J7TbmVJxqDW5xOUYv+xd/07aQcPnEBiVGLwupoeR3LMEAb0HwOAR/k4ePSdVutazFCTw2Z8BZ2pcbPxiE72dL6hHYSm0xwoqOBAQWPKyYvHJzWR8W6OtzSLveXHjILFyoT0S8JtouYcIiJMTF8KptjintKD+MpyWq0fH2Vn4ZjGHC67cso5WlTVan3NuUc7CE2nKK9x8+HBQn85Y0gco5Nj22yTdeQdanxG0qDo+OGMHDA+rDZqzj1jkjKJiB0CQIXPRXZW66MIgAlD4hgVkFlw7f7T/pzlmu5HOwhNh/H6FO/sLcDlMWLGCVF2Fo1tO5ufqipiV+FOoyDCxLSLsVqs4TZVc46xW+1kBOgz7T61vdUlr2CMOpZkJBPrNEaetS4v7+49pQX9zhO0g9B0mE+OlVBgrjqxiHDZxMFE2Nr+si84tpZiby0A1qhEMoZeGHY7Nd3DpNQFSKSxBybfXUnRsZYqr4E47VYumTCoIUEdJ0tr2BGwKk7TfYTVQYjIpSJySESyROSBIOfHichWEakXke91pK3m/7d338Fx3mdix7/Puw1YlEUvBMACEiIJkqIIUiQt6SxblmRRlERJJ5/l+Oyck4nGjjW55JJJ7DiTTCbJ3M1cyiVzntNpLorHPsuyzrJ0PPVuWVZhAXsHWNCJvqjbf/njXS4Acgmi7HJB4vnMYMh3d993H1HYfd5fe36Z0dI3xp7zU1dLX9o57KpCYxxu+zhxuLrmLrKd2ekKUWVYrjuXlVUTU14Pt/wGItOvmK4p8k7ZO+J3TX10DI6nLUY1M2lLECLiAH4M7ADqgW+ISP1lL+sH/gXw3+dwrrrORoIR3jjayaXWf02RN3khvsv4L/yWc6F43R1PLrfqntM3vY21XwWXfRNwJtDDWOvn1zgDttcWUxm/2YgZw+tHOgmEk6+lUNdHOlsQW4Gm+PahIeAFYNfkFxhjuo0xe4HwbM9V11csZnjzaBdjIfsDm+NxsGP9lXs8XCES4vDZtxN9yjUVDRRlF09/jrrhledUUFG2HrC/7A83v3bVhXOXOCxhx4ZKPC77a2k4EOGtYzoekUnpTBBVwOQ5bm3xx1J6rog8JSL7RGRfT4/Oo06Xz8/109o/BoAIPLCukpxpprReMtb6GcfH4wvpnFnctnJHOsNUC8jGlTshvkr+6Egbwfb91zzHl+3i/vqJwo1ne0ZpbNGqr5mSzgSR7NZyprcCMz7XGPOsMWaLMWZLaen0M2nU3JzvHeXzcxMbAG1dUcTSYu+1T4xGONT8WqJqa2nZeqrzl6UrTLXArCiqo6B4NWBXeT165tWke0VcblVZLpuWTmw9+/GZXtp1PCIj0pkg2oCaScfVQMd1OFelkH88zBtHu6aMO2xfMbMuomDHfo6NtNsHDjeb6x7Wqq2LiCUWDXUPg2W3NA8PnyN88eiMzv29utLE5IeYMbx2uEPXR2RAOhPEXqBORFaIiBt4Eth9Hc5VKRKJxnj1cEdioDAvyzmzcQeAWIwjZ14lZOxzC0vWsKKwLp3hqgWormQduUUrARiPRThxejfMYEzBYQk7b60k221Pnx4NRnn9SKeWBr/O0pYgjDER4GngLeAE8KIx5piIfFdEvgsgIhUi0gb8CfAfRKRNRPKvdm66YlVXMsbw/snuxO5wDkt4cMPMxh0Awl2HOTx83j5wOGnQ1sOi5LAcbFq1E+KLIg8Mnibae2pG5+ZnuXhwfWVifUT7wLjWa7rOZvZpnyNjzOvA65c99sykv3dhdx/N6Fx1/RxsHeRYx1Di+Iu3lLKkYIZrF2JRjp16hUDM7hLIK1zFKi2rsWitKb+N/QXLGOs/y2gszKkTL1Nf8gOYwQ3D0mIvd6ws4XdNvQAcaBmkNM/DuiVXrxisUkdXUqsrtPSN8dHp3sTx2sp8Nk5TwvtyoY5GDgw12weWk011D2lZjUXMZbnYuPJBsOyvm30DJ4lcnHmHwO3LC1k5qV7Teye6dRHddaIJQk0xOBbitSOdia1DK3xZ3Lu2bObdQ9EIh06+zHi89ZBbXMeaso3pClfdINZXbiW7cAUAI7EQx06+NKMZTWDXa/rqunJKcu29JqIxw6uHOxgOXL58SqWaJgiVEAhH2X1oYlA61+Pk4Y1LcDpm/msy3voph0Za7AOHm9tveRSnldaeTHUDcDlcbK7blZjR1OhvItzROOPzPU4Hj2ysmjJovftQR6JgpEoPTRAKuHRX1pnYV9ppCQ9vXDLt/g5XiAQ5cHp3YuZSQWk9q0vXpyNcdQOqL99EbrE9k208FuHQqZchOvOpqz6vi50bKrHirdnuoSBvHO0kpjOb0kYThMIYw7snLiZWSgPct6782kX4LjN67iOOjMWXqziz2HbLo1iiv2LK5rSc3H7Lo4nV1QeHLxBo/WxW16gp8vLlNRMLYs/2jPIbndmUNvrpVew518/xSTOW7lhZzJqK/NldJOBnb/M/TKyarthIbXz7SaUuWV26noJSu+5myERpPPUyhEZndY1bqwvYsnyiSOTBlkEaW7Q8eDpogljkjrb7+aR5oozGuiX5bF1RNM0ZyfWc/AdOBOIzn9w5bKt7VNc9qCtYYrF19WOJSq+HxzsZPP3mrK9z16oS6sonZjZ9dLqHU13D05yh5kITxCLW1D3CuycuJo6XFnn5ytryWX+xm4EWPm79cKJi69K7qPFpzSWV3MrCW6is2grYZTQ+Pv82DF+8xllT2TObKlhSYHeDGgNvHeviQt/sWiNqepogFqnW/jHeODKxt0NZvoeHNlbimEkZjcmMoenYC3RG7A+m5JZyV90ubT2oqxIR7rrlccRrt1RbwkNcOPbijEpwTOZyWDyysYriKdNfO+mK73ao5k8TxCLU5Q+w+1AHkfjsj0Kvi8c2VV1z29Bkwu37+bQvvuhJLG5dtZPCrGtvIqQWt9KcMtbU3p9YTf1x936i3cdnfZ1st4NHN1Ul9rQORWK8fKCdnuFgSuNdrDRBLDLdQwF+faAtMX881+PksYZqvO45rFUIjXLgxIuMxOypsdlFtWxZ+qUURqtuZtuW34fbtxQAfzTI4aPPQ3j2d//5WS4eb6hOrJEIhKP8urGNvhFNEvOlCWIR6RkO8lJjO8GwnRyy3Q4ea6jCl+2a0/X6jr9M42ibfeDMYtuaJ/A4PKkKV93kvC4vW1Y/Dg67i2jv8Hn8p1+b07WKctw8tqkKt9P+ShsLRXmpsY2B0en3wlbT0wSxSHQPB3ipsS2xSjrL5eDxhipKcuf2hR7rOcUHrR9MlOSo2c7a0ltTFq9aHDZUbKa4sgGAiInx4dk3Mf3n5nSt8vwsHm+YSBKjwSi/2t9GvyaJOdMEsQh0+QO8tL+d8fh+0h6XxeMNVZTlzW4hXEI4wOHDP6U7Yi+ss/Iq+NKar+nAtJo1h+Xgy/VPIjklALRHRjh++GcQnVudpUpfNrtuW4LLYf8ujgQj/N2+Vh2TmCNNEDe5toGxKS0Hj8visU1VlOfPMTkA/lOvsiex14OLLWufoChr9msnlAIoyyln4+pHE3WaPh08zciZt+Z8vepCL7tum9rd9Kv9bVwc0tlNs6UJ4iZ2rneUVw60Jwaks90OnmioptI3w30dkoh1n+D9c28Sia+YLlmymU2V21MSr1q8bq/5Ir7yDYC9wvqDpt2YvuY5X6+myMujk8YkAmE7SUwuJ6OuLa0JQkQeEJFTItIkIj9I8ryIyP+JP39YRBomPXdeRI6IyEER2ZfOOG9GR9v97D7YQThqjxHkeBw8sbmasnm0HAj42X/wuSlrHr609knd60HNm8ty8eV134Rse4p0a3iYQweem3UZjsmqCrJ5YnM1WS779/PSFFhdcT1zaUsQIuIAfgzsAOqBb4hI/WUv2wHUxX+eAv7qsue/bIy5zRizJV1x3myMMXx2to93jl9MDCDnZ7v42uaaOQ9IAxCL0X7ob9k31m4fO7O4vf5JynLKUxC1UrAkr4rb1j6RmNX02cg5uo+8MOsFdJOV52fxxObqRFXiaMzw+pFOrd00Q+lsQWwFmowxZ40xIeAFYNdlr9kF/NTYPgMKRKQyjTHd1CLRGG8du8ink2orleZ5+PrtNRTmuOd17fGzH/Bu1+d2OQ0Rliz7Ig1LtGtJpda2mrspq/kCYJfheKf9t4RaPpnXNUvzPPzB7TUUTfoM/OZUD++fvEhUS4VPK50JogponXTcFn9spq8xwNsisl9Enrram4jIUyKyT0T29fQs3rK/o8EILzW2caJzoirr0iIvX9tSPbs9HZKIdZ/kvZO/ZDRmzyzJKlrFvWu+pqW8Vco5LAf31X8Dd6Fdy8sfDfLB0Z9j+s7O67q+bBdfv70mUbsJ4FCrn5cPtCcmcKgrpfMTnmzO4+XperrX3GmMacDuhvq+iHwx2ZsYY541xmwxxmwpLS1N9pKbXpc/wC/2tNAxODFLY92SfB6dY/mMKUa6+azxr2kJx/ttswu459Y/ItedO/15Ss2Rz+Pj7g3fAU8eAM2hAfY1PgNj/fO6rr32p5rVFXmJx1r7x/jFnhadBnsV6UwQbUDNpONqoGOmrzHGXPqzG3gZu8tKTWKM4VDrIC/ua2U4YO/MJQJ3ry7lvvry2Rfeu1xojJP7/pqD4132sSuLTfVfZ3lB7TwjV2p6dcWrWb/2cXDGV1mPtNK8/1mIzO+L3OWw2LG+gjtWFiceGxwL88u9LRzr8M/r2jejdCaIvUCdiKwQETfwJLD7stfsBr4dn820HfAbYzpFJEdE8gBEJAe4HziaxlhvOMFIlDePdvH+ye5EP6rHZbHrtioalhbOf9FaNEzngZ/w4eAJ+9hysHzlV9lec/c8I1dqZu5cdh/VK+6BeFfme31H6Dn4M4jNr0tIRNhWW8zDGysT02DDUcPbxy7yzvGLhKO6z/UladtN3hgTEZGngbcAB/CcMeaYiHw3/vwzwOvAg0ATMAZ8J356OfBy/EvOCTxvjJn9riI3qY7Bcd482oV/fGK1aWmeh4duraTAO7/BaABiUfoO/pzXuz5NzIQqrtrKvbc8rqul1XXjsBzcv+brvDTag79jPxET47X2j3jU5aVgw5Ngze/+dlVZHoVeN68dmdiL/Wi7n07/OA+sq5jflPCbhJh5TCFbaLZs2WL27bt5l0xEY4Y95/rZc64/8cUNsKHKx92rS3E5UtAgjMUYPPpLXjn3BmPxQens0jX8/tZ/Rb57ltuQKpUCA4EBXvr8fxDqtxfO5VpuHlu1i7z6xxLlwucjFInx3omLnJy0PsJhCV9YWczmpYVY8+2qXeBEZP/VlhLoNJQbxMUheyD6s7N9ieTgcVns2FDBvfXlqUkOxjB04hV2n3szkRzcRSvY2fA9TQ4qYwqzCnmw4Xs4ffZw5UgsxO7m3YydfmNeayQucTstHlhfwX315Ykup2jM8PGZXl7c10rvIi4brgligQtHY/z2TM8VMy2qCrP55rZlrKlI0Rd3LIr/yAvsbvr7xP4OTl81Oxuepsxblpr3UGqOluRV8cDm72Pl2Qsz/dEgr5x6keHjr0Bs/mMGIsL6Kh//aOtSKnwTXUud/gDPf97CJ829RBbh2IR2MS1QxhjOdI/w0emexAwlAJfDbvpuqklh0zcapufgz3i17QPGY/Z7WXnl7Nz6J9TkL03NeyiVAmcHmnlr719gRnsByLFcPLzsfopu/QakqORLNGbYe97uyp28kM6X7eLu1aXUluTcVGNx03UxaYJYgLqHA3x0uveKwmI1RV7uXVuWmoHoS4IjtB14jje79hAy9uwQh6+arzZ8T6ezqgWpqf8U7zY+Q2zYnn7tsZzsrLyLik3fBtfcC1Ferm8kyLsnLk5ZXwSwvMTLF+tKKZ5P6ZoFRBPEDcI/FuaT5t4pg2UAXreDO1eVsG5JfkrvXIy/nSP7n+F3Q81c+j1wF9Wys+H7VOYtSdn7KJVqrUOtvLn/LwkPtgDgEIvfK1hN/ZbvQm7qukSNMRxu8/NJc9+UFdciUF+Zz/aVxeRnzW1HxoVCE8QC5x8Ls+9CP8c6hqY0aS0RNtb42F5bnKhImRLGEO5o5MND/48zwXjdJhFyStfyUMM/pzi7ePrzlVoAukcv8ur+HxPoO5N4rD67lLs2/jOclRtS+l7joSifNPdypN0/ZVzcaQnrq31sXlZ4wyYKTRALVO9IkH3nBzjVNTxl2ipAbWkOd64qmV8F1mRCo3QffZH32j5iIBpvOltOSmu+wAPrv0WeO2/685VaQPxBP28c+Qn97XsTC+hKHNncs+xeSuofS2mXE9jdv79r6uV879TuX4cl1Ffms3lZ4bwLY15vmiAWEGMM5/vGONAywIW+KzcvqSrI5s66EqoKUvuLDRC5eIL9h39C42hroksJt5f61Y9yV+0OnFba1k0qlTbhaJgPm3Zz5sxrEB4H7Nb3lrzlNGz8DlZJXcrfs7V/jI+beunyX7lLXW1pDptqCqkpyr4hBrM1QSwAI8EIxzuGONrun7IC+pKaIi/bVhRRXZiGX6rRPlqP/x0fd34+0WrAnsZ61/o/pL5sY2rfT6nrzBjDka79fHrs50Tjg9dgtybuqrqTJfWPJzYjSuV7XugbY8+5ftoHx694vtDrYn2Vj7WV+eTMs6JyOmmCyJBgJMrZnlFOdQ1zvm/0ijU9IrCyNJfNywpZkoYWA6FR+pvf4ZPmN2kJTypE5nRTWXMn96z9Oj6PL/Xvq1SG9Af6ef/Y83S374HoxI1YraeQL6zciW/Fl8HtTel7GmNoGxinsWWAsz1X7oBnibC8xMuainxqS3NSs6g1hTRBXEeBcJRzvaM094xwvnc0seXnZB6XxbolPm6rKcCXnYaBreAIvc1v03jhfZqD/RPdSSK48qvYtvYP2FCx5YZo/io1WzET42D7p+w9+asprQkR4RZPMZtW3E9R7T3gzkn5ew+MhjjYOsjxzqHEXvCTuZ0WK0pyWFmay/IS7/zL8aeAJog0MsbQMxzkQv8Y53tH6RgMXDHgfElNkZf1VfmsLM1N/V2EMcQGWzh/9h2Ode2jNTS1dLFk+Vhdez/bar9Kjiv1HwylFprh0DCfNb3OmfPvQXDq1PHl7kLWLdlKTe29WL7qlL93KBLjTPcwx9qHknY/gT2wXVWQzfISL0uLcijJdWfkpk0TRApFY4a+kSBtg+O0D4zTPjjOeOjq5YdL8jysqcjjlvK8tLQWzFg/fW2f09T+KaeGzid2fUvw5FKzZCvb63ZRmqMlM9Ti0zXSyeenX6G9cx+Epk4MybXcrCmoZVXVHRRV3Q7ZBSl//8GxECe7hjnVNUz/aOiqr8vxOKgq8FJVmM2SgixKcjzXpVCgJog5isUM/WMheoaDdA8HuegP0D0cSNptNFmFL4uVpbmsLM1J/WrLWIyIv5Wuzv20dh/m7NB5/NEri4lJVj611XfQULtDE4NS2Imisek1znd8DsGRK54vdGRR66ulpvxWyisbcORXp6Ra7CXGGHpHQjT3jNDcM0L30PRFAN1Oi7I8D+X5WZTleyjN9VDodac8aWiCuIZAOIp/PMzAWIjBsTD9oyH6RkMMjIZmtKl5ttvBsiIvS4u9LCvOmfce0FOEAwT9rXT3Hqer7xRd/vN0hvxETJLCYZaFN7+GNUu/yNrqO3UAWqkkBsb7OdH2MSdbf0tgqAOSfJZc4qDSU0CFbwXlRXWUl67D7asGZ+pu+IYDYS70jXGhb4yW/rEZ7Y3ttITCHDfFOW6KctwUeN0Uel34vK45j2dogpjGC3ta6Ewyl3k6eVlOKn3ZVBVmU12YTXHOPPsOjcEERxgf7WZoqIXB4XYGRzrpH+mkb7yP4djVm6WIhctbzPKyjaysuYtlhatwpKhomVI3s0gswoX+0zS1/o4L3YeJjPdNWz4833JTnF1CUW4lBbmV+PKqyPctJTunDHHnzKu1YYyhZySY6Lbu8gemFOmciaVFXn5/8+zHU6ZLEGmdnCsiDwD/G3tHub8xxvzZZc9L/PkHsXeU+yNjTONMzk0Vj2v6weK8LCeleR5K8zyU5WVR4cuavoVgDCYaJhoZJxwaJRQeJRwaIRgaIRj0EwiNEAgNMR7wMxYaYiw0zEhwiJFokGiyVkEyTje+vCpqStZRU7mZmsJVushNqVlyWk5WltSzsqSecCxMa99pWroaaes9xtBwx5RpsgBDsRBDox2cG+2Ai5OuIxY5Dg95Hh9edz7Znjy8nnyy3D48nlw87nw87lxc7lzc7lxcLi8OZzbicCWSiohQlpdFWV4Wm5ba6zWGA2EuDgXoHgrSMxKkeyjISPDqSSOl5Xgu/bel/IpxIuIAfgzcB7QBe0VktzHm+KSX7QDq4j/bgL8Cts3w3JQID7xFpPccHpfgcTnIcghZLsHjtP+0MMSGYvj9MQZMlBMmRiwWJWYiRGNRYrEokViEqIkQjYWJxCKETYyUtcxEsNy5FOVUUFG4ivKStVQW30K+dh8plTIuy0Vt6TpqS9cB4A8M0tl3kou9J7g4eJa+kU5MeCxpCyNiYvgj4/gj4zDadcXzyYgILrFwWk4clguHOHFaTizLgcNyYokDy3JgiQOHOPCIsEwsYh4hEDYEI4bxsCEYiREMxwhEDBH3auBbqfxnSWsLYivQZIw5CyAiLwC7gMlf8ruAnxr72/QzESkQkUpg+QzOTYkszuHNakaAaAxGYzB65ULn9LOcuN1e8rOK8HnL8OVWUORbTlFBLYXeEu02Uuo68mUV4Kvazpqq7YDdHTU41kvfQBMD/hYGR7sYGuvGP95PODyWqAM1U8YYQiZKKBYF5rdjndMJOU5wS+rLj6czQVQBrZOO27BbCdd6TdUMzwVARJ4CngJYunT2m9s4LIuUTyQTC8ty4nS48Dg8uJ1ZuB1ZeFxePO4cst15ZHt8ZGcV4M0uISe3nLzsYtyOG6vIl1KLhdNyUpJbQUluBdRMfS4YDTI82svo6EXGA/2MjQ8wHvQTCI8QDI0SDI8RigYIRgKEokEi0TCxWCTp4PhcCeC0Ur9CO50JItn37uXts6u9Zibn2g8a8yzwLNiD1LMJEKB++VdYWrYBEQtBwLKwsOxjsZB4M88Sh/13y5n4cVhuHA43lsOF0+nB6cyyfxweveNXapHwODx48qsoya+a8TnRWJRINEgkPG7/GQkQi0aIRkNEo0FiJkosFiEWi2BiUWImSjSeVEz8J0YMYvFjDL407OGSzgTRxtRcWw10zPA17hmcmxIrlt+djssqpdRVOSwHDsuLx5XaulCpls6qUXuBOhFZISJu4Elg92Wv2Q18W2zbAb8xpnOG5yqllEqjtLUgjDEREXkaeAt7qupzxphjIvLd+PPPAK9jT3Ftwp7m+p3pzk1XrEoppa606BfKKaXUYjbdQrmFVZhcKaXUgqEJQimlVFKaIJRSSiWlCUIppVRSN9UgtYj0ABfmeHoJ0JvCcFJF45odjWt2NK7ZuRnjWmaMKU32xE2VIOZDRPZdbSQ/kzSu2dG4Zkfjmp3FFpd2MSmllEpKE4RSSqmkNEFMeDbTAVyFxjU7GtfsaFyzs6ji0jEIpZRSSWkLQimlVFKaIJRSSiWlCSIJEfk3ImJEpCTTsQCIyH8RkcMiclBE3haR1O8MMgci8ucicjIe28siUpDpmABE5GsickxEYiKS0SmJIvKAiJwSkSYR+UEmY5lMRJ4TkW4ROZrpWCYTkRoR+UBETsT/H/5xpmMCEJEsEdkjIoficf3nTMd0iYg4ROSAiLya6mtrgriMiNQA9wEtmY5lkj83xtxqjLkNeBX4j5kOKO4dYL0x5lbgNPDDDMdzyVHgceCjTAYhIg7gx8AOoB74hojUZzKmSX4CPJDpIJKIAP/aGLMW2A58f4H8mwWBe4wxG4HbgAfie9gsBH8MnEjHhTVBXOl/Af+Wq2xxmgnGmKFJhzkskNiMMW8bYyLxw8+wd/7LOGPMCWPMqUzHAWwFmowxZ40xIeAFYFeGYwLAGPMR0J/pOC5njOk0xjTG/z6M/cU3870808TYRuKHrvhPxj+HIlIN7AT+Jh3X1wQxiYg8ArQbYw5lOpbLich/E5FW4JssnBbEZP8EeCPTQSwwVUDrpOM2FsCX3Y1CRJYDm4DPMxuJLd6VcxDoBt4xxiyEuP4C+4Y2lo6Lp3NP6gVJRN4FKpI89SPg3wP3X9+IbNPFZYz5e2PMj4AficgPgaeB/7QQ4oq/5kfYXQM/vx4xzTSuBUCSPJbxu84bgYjkAi8B//KyFnTGGGOiwG3xsbaXRWS9MSZjYzgi8hDQbYzZLyJfSsd7LLoEYYy5N9njIrIBWAEcEhGwu0saRWSrMaYrU3El8TzwGtcpQVwrLhH5x8BDwFfMdVxUM4t/r0xqA2omHVcDHRmK5YYhIi7s5PBzY8yvMx3P5YwxgyLyIfYYTiYH+e8EHhGRB4EsIF9E/tYY84epegPtYoozxhwxxpQZY5YbY5Zjf7gbrkdyuBYRqZt0+AhwMlOxTCYiDwD/DnjEGDOW6XgWoL1AnYisEBE38CSwO8MxLWhi3539X+CEMeZ/ZjqeS0Sk9NIsPRHJBu4lw59DY8wPjTHV8e+rJ4EzHvJeAAABaElEQVT3U5kcQBPEjeLPROSoiBzG7gJbEFP/gL8E8oB34lNwn8l0QAAi8piItAFfAF4TkbcyEUd8AP9p4C3swdYXjTHHMhHL5UTkF8CnwGoRaRORf5rpmOLuBL4F3BP/nToYv0POtErgg/hncC/2GETKp5UuNFpqQymlVFLaglBKKZWUJgillFJJaYJQSimVlCYIpZRSSWmCUEoplZQmCKWUUklpglBKKZWUJgil0kREbo/vlZElIjnxfQTWZzoupWZKF8oplUYi8l+x6+RkA23GmD/NcEhKzZgmCKXSKF6DaS8QAO6IVwRV6oagXUxKpVcRkItdsyorw7EoNSvaglAqjURkN/ZOciuASmPM0xkOSakZW3T7QSh1vYjIt4GIMeb5+P7Un4jIPcaY9zMdm1IzoS0IpZRSSekYhFJKqaQ0QSillEpKE4RSSqmkNEEopZRKShOEUkqppDRBKKWUSkoThFJKqaT+Pw9PUj/XT5KxAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(support, init_untwist_dist, lw=3, alpha=.5, label = f\"t=0 untwisted distribution\")\n", + "plt.plot(support, init_twist_dist, lw=3, alpha=.5, label=f\"t=0 twisted distribution\")\n", + "plt.plot(support, final_dist, lw=3, alpha=.5, label=f\"t=9 (last) final target distribution\")\n", + "plt.xlabel(f\"x\")\n", + "plt.ylabel(f\"density\")\n", + "plt.title(f\"initial and final target distributions\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "so it would appear that we can recover the exact posterior since the prior family includes the posterior.\n", + "\n", + "In a physics context, this is akin to annealing between an oscillator of spring constant $K_0$ to a spring constant of $2K_0$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, I'll attempt to anneal between two springs of a constant $K_0$ but from $r_0=0$ to $r_1=2$" + ] + }, + { + "cell_type": "code", + "execution_count": 99, + "metadata": {}, + "outputs": [], + "source": [ + "def new_potential(pos, parameter):\n", + " mod_parameter = parameter[0]\n", + " dim = pos.shape[0]\n", + " r0, r1 = jnp.array([0.0]), jnp.array([2.0])\n", + " k0 = jnp.diag(jnp.ones(dim))\n", + " delta0, delta1 = pos-r0, pos-r1\n", + " return 0.5*((1. - mod_parameter) * jnp.dot(jnp.dot(delta0, k0), delta0) + mod_parameter*jnp.dot(jnp.dot(delta1, k0), delta1))\n", + "\n", + "#make a dummy potential for fun\n", + "def dummy_potential(pos, parameter):\n", + " return 0." + ] + }, + { + "cell_type": "code", + "execution_count": 170, + "metadata": {}, + "outputs": [], + "source": [ + "#wrap an initializer\n", + "def get_test_potential_initializer(dimension):\n", + " potential = new_potential\n", + " mu, cov = jnp.zeros(dimension), jnp.diag(jnp.ones(dimension))\n", + " dG = 0.\n", + " return potential, (mu, cov), dG\n", + "\n", + "#define some _vanilla_ A, b fns\n", + "A_fn = lambda x, params: params\n", + "b_fn = lambda x, params: params\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "next, I'll anneal between a gaussian with a covariance of 1 centered at 0 and anneal to the same gaussian, except centered at x=2 (again, over 10 steps)" + ] + }, + { + "cell_type": "code", + "execution_count": 268, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "DEBUG:csmc:instantiating super\n", + "DEBUG:csmc:define A, b functions\n", + "DEBUG:csmc:getting the parameter lengths\n", + "DEBUG:csmc:getting x dimension\n", + "DEBUG:csmc:getting the terminal time\n", + "DEBUG:csmc:getting the maximum number of twists\n", + "DEBUG:csmc:getting prop_params\n", + "DEBUG:csmc:getting init_params\n", + "DEBUG:csmc:getting model_params\n", + "DEBUG:csmc:getting y\n", + "DEBUG:csmc:define twisting parameters and parameter caches\n", + "DEBUG:csmc:define psi twisting functions\n" + ] + } + ], + "source": [ + "N,T,Dx = 10000, 10, 1\n", + "potential, (mu, cov), dG = get_test_potential_initializer(Dx)\n", + "cov = cov[0]\n", + "potential_params = jnp.linspace(0,1,T)[..., jnp.newaxis] #define the potential parameters\n", + "\n", + "prop_params = {\n", + " 'potential_params': potential_params,\n", + " 'forward_potential_params': potential_params,\n", + " 'backward_potential_params': potential_params[1:],\n", + " 'dt': dt\n", + " } #potential parameters dictionary goes here\n", + "init_params = {'mus': jnp.array([mu]), \n", + " 'covs': jnp.array([cov]),\n", + " 'mixture_weights': jnp.array([1.])} #define the initial params\n", + "model_params=None #there are no model parameters\n", + "y=jnp.arange(T) #the y data aren't important. just know that there are the same number of them as there are data\n", + "rs = random.PRNGKey(6) #define a random number \n", + "\n", + "#create the csmc container object\n", + "smc_obj= StaticULAControlledSMC(N, \n", + " potential=potential,\n", + " forward_potential=potential, \n", + " backward_potential=potential, \n", + " A_fn = A_fn,\n", + " b_fn = b_fn,\n", + " params_len=1,\n", + " Dx=1,\n", + " T=T,\n", + " max_twists=10,\n", + " prop_params=prop_params,\n", + " init_params=init_params,\n", + " model_params=model_params,\n", + " y=None\n", + " ) #create a controlled smc class" + ] + }, + { + "cell_type": "code", + "execution_count": 269, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "DEBUG:csmc:generating SMC variable initializer function\n", + "DEBUG:csmc:generating simulation propagation function\n", + " 0%| | 0/9 [00:00" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "sns.distplot(controller.bootstrap_logZs[0], label=f\"uncontrolled SMC\")\n", + "sns.distplot(controller.bootstrap_logZs[1], label=f\"twist 1\")\n", + "sns.distplot(controller.bootstrap_logZs[2], label=f\"twist 2\")\n", + "sns.distplot(controller.bootstrap_logZs[3], label = f\"twist 3\")\n", + "#sns.distplot(controller.bootstrap_logZs[4], label = f\"twist 4\")\n", + "#sns.distplot(controller.bootstrap_logZs[5], label = f\"twist 5\")\n", + "#sns.distplot(controller.bootstrap_logZs[6], label = f\"twist 5\")\n", + "plt.vlines(dG, 0,1000, label=f\"true logZ\")\n", + "plt.ylim(0, 100)\n", + "plt.legend()\n", + "plt.ylabel('counts')\n", + "plt.xlabel(f\"logZ\")\n", + "plt.title(f\"bootstrapped logZs\")" + ] + }, + { + "cell_type": "code", + "execution_count": 274, + "metadata": {}, + "outputs": [], + "source": [ + "support = np.linspace(-4,4,100)[..., np.newaxis]\n", + "init_untwist_dist = np.exp(multivariate_normal.logpdf(support, mean = mu, cov = cov))\n", + "_, _, (twist_mu0, twist_cov0) = get_twisted_gmm(controller.csmc_object.init_params['mixture_weights'], \n", + " controller.csmc_object.init_params['mus'],\n", + " controller.csmc_object.init_params['covs'],\n", + " controller.csmc_object.A0,\n", + " controller.csmc_object.b0)\n", + "init_twist_dist = np.exp(multivariate_normal.logpdf(support, mean = twist_mu0[0], cov = twist_cov0[0]))\n", + "final_dist = np.exp(multivariate_normal.logpdf(support, mean = 2., cov = np.array([1.])))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 276, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 276, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3gc1dXwf2d3tepdcpVtybLci9wxlm3MCwaCDQ4ltIAJ4SUmgff9yBdKQkJJPoJJnJAQQoCQxKQQSBwgpiVAwNimuYCNe5dtWZZVrF633O+P2V3NrlbSrqzVrlbzex492il37pmZ3Tlzzz1FlFIYGBgYGBj4Ygq3AAYGBgYGkYmhIAwMDAwM/GIoCAMDAwMDvxgKwsDAwMDAL4aCMDAwMDDwi6EgDAwMDAz8YiiIKEFEdovIeb2xr4i8JSIrAjxWsYhcEJiUPUdEzhORki62zxeRgyLSICLLgzmHs+m3P6K/ZyLyPRF5rheP3SAio12f14jI/+vFYz8tIj/oreMZdI8l3AIY9A5KqUk92VdEHgLGKKW+qtt+Se9K1yf8EHhSKfVL1/KrfdGpiCigQCl1qC/689P/euDPSqkePeSVUj/uzX6UUkk9kcNPfzcDtyqlinTHXtkbxzYIHGMEYRAtjAJ2h1uIYBCNqPgNiojxshmFRMWX06CD2eAhEfmbiPxRROpdJqVZvvuKyMXA94BrXKaBHa7t60XkVtfnfBF5T0SqRKRSRP4iImkBynSpiHwuInUicsI1WnFvyxURJSIrROS469j367bHu0wU1SKyB5jdRT+HgdHAa67ziPU5h5tFZJOIrHYd76iIXKJr/zUR2eu6VkdE5BsBnt8G18cdrn6vEZF0EXldRCpcfb0uIjm6NutF5BER+RBoAkaLSJ6IbHD1/66I/FpE/qxrc46IfCQiNSKyw20eFJFHgAXAk67+n+xEzhtF5JjrHt7vs+0hd18iEicif3btVyMiW0RkcGf9uO7ft0TkIHBQt26MrossEXnHdW4fiMgo137u+2/RybJeRG4VkQnA08A8V381ru1eJisR+W8ROSQiZ0RknYgM021TIrJSNLNjteuaimvbGJcsta7v3Uvd3+0BilLK+IuCP6AYuMD1+SGgBfgSYAYeBT7pYt8/+xxrPdrwHmAMcCEQC2QDG4Bf+DuWH5nOA6agvYhMBU4Dy13bcgEF/BaIB6YBrcAE1/ZVwEYgAxgB7AJKAjl/P+dwM2AD/tt1PW4HSgFxbb8UyAcEWIT24J6hO4eu+lVoJjr3ciZwJZAAJAN/B171kes4MAnNxBsDfAysBqxAEVDnvifAcKDKdS9NrntRBWT7nmcn8k0EGoCFrnv4c8Du7/4D3wBec8luBmYCKZ314zr3d1z3KN73egBrgHpd378ENvncf0sX92yTT39rgP/n+nw+UAnMcB37V8AGH9leB9KAkUAFcLFr21+B+13XMw4oCvfvN1L/jBFE9LJJKfWmUsoB/AntARw0SqlDSql3lFKtSqkKtAfMogDbrldK7VRKOZVSX6D9MH3bPqyUalZK7QB26OT8CvCIUuqMUuoE8ERP5NdxTCn1W9f1eB4YCgx2yfmGUuqw0vgAeBvtjTlolFJVSql/KKWalFL1wCN0POc1SqndSim7S47ZwANKqTal1CZgnW7frwJvuu6lUyn1DrAVTWEEwlXA60qpDUqpVuAHgLOTfW1oCm6MUsqhlNqmlKrr5viPuu5Rcyfb39D1fT/aqGBEgLJ3xQ3A75VSn7mO/V3XsXN1+6xSStUopY4D7wOFrvU2NJPkMKVUi+uaG/jBUBDRS5nucxMQJz2wE4vIIBF5UUROikgd8GcgK8C2c0XkfZe5pRZY6aetr5zuSc5hwAndtmPByt5ZP0qpJtfHJJecl4jIJy5TRQ3awzegc/RFRBJE5BmXSacObcSVJiJm3W768xoGnNHJ5Lt9FHC1y+RT45KvCE2xBILXdVRKNaKNQPzxJ+DfwIsiUioiPxGRmG6OfyLQ7UqpBuCMS6azZRi674Tr2FVoIy43nX237kEbLW4Wzfx6Sy/IE5UYCsKgu3S+j7r2maqUSkF7o5UAj/0C2tvwCKVUKppdOdC2p9BMS25GBtguKEQkFvgHmolnsFIqDXiTwOX05f8C44C5ruu10N2Vbh/9NT8FZIhIgm6d/rxPAH9SSqXp/hKVUqv8HMsfXtfR1U+mvx2VUjal1MNKqYnAucBS4KZu+umuf33fSWjmqFKg0bVaf95DgjhuKZrydB87Ee28TnbTDqVUmVLqv5VSw9DMak/5zJsYuDAUhMFpIFc696ZJRrNh14jIcODuII6djPZ23CIic4Drg2j7N+C7rknfHODOINoGgxXNhl0B2F2T10uCaH8abYLcTTLQjHa9MoAHu2qslDqGZjJ6SESsIjIPWKbb5c/AMhG5SETMronk83QT3779+7IWWCoiRSJiRXMH9nuvRWSxiExxjXbq0EwxjgD76Ywv6fr+EfCpUuqEy1x5Eviq67xuQZsHcnMayHG188cLwNdEpNCl5H/sOnZxdwKJyNW661eNpowcXTQZsBgKwuDvrv9VIvKZn+0Po00E1gJvAC8HcexvAj8UkXrgAbSHfqA8jGZCOIo2J/CnINoGjGue4H/QZKtGU2LrumzkzUPA8y7zz1eAX6BNulcCnwD/CuAYNwDz0Ewk/w94CW3CHtf8y+Vo3mYVaCOKu2n/7f4SuMrlqdNhnkYptRv4FtoD9ZTrHDsL/BuCplDqgL3AB2gKqtt+uuAFNCV5Bm3S+wbdtv92nUsV2qT9R7pt76G5LZeJSKWf8/oP2nzKP1znlQ9cG6BMs4FPRaQB7V7/r1LqaBDnNGBwe3EYGBhECC63y31KqS5HHwYGocYYQRgYhBkRmS1avIlJtNiUy+mjSHADg64woh8NDMLPEDTTXSaa+ed2pdTn4RXJwMAwMRkYGBgYdIJhYjIwMDAw8EtUmZiysrJUbm5uuMUwMDAw6Dds27atUimV7W9bVCmI3Nxctm7dGm4xDAwMDPoNItJplgLDxGRgYGBg4BdDQRgYGBgY+MVQEAYGBgYGfomqOQh/2Gw2SkpKaGlpCbcoBgYRSVxcHDk5OcTEdJe41WCgEfUKoqSkhOTkZHJzc3EVlDIwMHChlKKqqoqSkhLy8vLCLY5BhBFSE5OIXCwi+11lAe/rYr/ZIuIQkauCbdsdLS0tZGZmGsrBwMAPIkJmZqYxwjbwS8gUhCtl8K+BS9DKHl4nIhM72e8xtEIlQbUNQpaeNjUw6DFKKexOJ3aH9udwOonEzAWR9PtosjVR1lhGTUsNNqct3OIMeEJpYpoDHFJKHQEQkRfRkpDt8dnvTrSUvbN70NbAIOJwOhXNNgfNbQ4cPgrBbBLiY8zEx5gxmSLnwRxOGtoaOFRziKO1RylrLEPpagXFW+IZlz6O6YOnE2+JD6OUA5NQmpiG412OsATvcoC4CtB8Ga3SWFBtdce4TUS2isjWioqKsxa6t6mpqeGpp54Kqs3Ro0eZO3cuBQUFXHPNNbS1tfWKLK+++ip79gSuY59++mn++Mc/drp9/fr1fPTRR51u74zc3FwqKzuk+PdizZo13HHHHb0ix7p161i1SivAdvPNN7N27dqg5P3xj3/stXzuuef63U8pRX2LjYqGVhpa7R2UA4DDqWhotVPZ0EpDiy0iRxR9hVKK3VW7eWHfC3xU+hGnGk95KQeAZnsz2yu285e9f2FL2RZsDmNU0ZeEUkH4ez3y/TX8ArjXVUg+2LbaSqWeVUrNUkrNys72Gy0eVnqiIO69917uuusuDh48SHp6Or/73e96RZZgFcTKlSu56aabOt3eUwURLGcjh91u57LLLuO++3o8jdVBQfjry+lUVDfZaGrz/ioLYBbBLNKh5mhjm4PaZhvOAagkWuwt/PvYv/ngxAfYnXbPekHIjMskKSbJy/TV5mhjS9kW/nHwHzTaGv0d0iAEhFJBlOBdWzcHrY6snlloBdKLgavQasMuD7Btv+C+++7j8OHDFBYWcvfd3VfrVErx3nvvcdVV2nz9ihUrePXVjqUB9G/YAEuXLmX9+vUAJCUlcf/99zNt2jTOOeccTp8+zUcffcS6deu4++67KSws5NNPP2XmzJkA7NixAxHh+PHjAOTn59PU1MRDDz3E6tWrAXjiiSeYOHEiU6dO5dprr6W4uJinn36axx9/nMLCQjZu3EhFRQVXXnkls2fPZvbs2Xz44YcAVFVVsWTJEqZPn843vvGNTt+a//CHPzB27FgWLVrkaQsELcfNN9/Mt7/9bRYvXsy9997b4Vq9++67LFiwgLFjx/L66693eT3vu+8+mpubKSws5IYbbvBcX/e9uvvuu5k8eTKTp0zh7397CYAPN27gikuXsPLmG1g4p5D/XXkLmUlWspJjSY6zYNGZllrtTs40tmF3OLv4VkQXdW11/P3A3zlSc8SzLj0uncUjFrNi0gquGX8NN026iZVTV3JR7kWkx6V79jvTcoaXD75MbWttOEQfcIRyDmILUCAieWi1Z6/FpyaxUsrjVycia4DXlVKvioilu7Y94fF3DpztITrlrgvH+l2/atUqdu3axfbt2wGor69nwYIFfvd94YUXGDRoEGlpaVgs2q3Jycnh5Mlu67B70djYyDnnnMMjjzzCPffcw29/+1u+//3vc9lll7F06VKP8mlpaaGuro6NGzcya9YsNm7cSFFREYMGDSIhIcHrmKtWreLo0aPExsZSU1NDWloaK1euJCkpie985zsAXH/99dx1110UFRVx/PhxLrroIvbu3cvDDz9MUVERDzzwAG+88QbPPvtsB5lPnTrFgw8+yLZt20hNTWXx4sVMnz7d7/XsTo7f/e53HDhwgHfffRez2cyaNWu8jlFcXMwHH3zA4cOHWbx4MYcOHer0Wq5atYonn3zSc//0vPzyy3z++ee8s+lTqioruXhxEfPmFxEfY2bXFzvYvXs3w4YNY/78+Xz44YcUFRWRYLUQH2OmodXuGW04nIrqpjYyEq2YTdEdu9pib+H1w69T31bvWTc5azLzhs0jxuQdhyEi5Kflk5eax96qvWw4uUEz47XV88rBV1iav5Ss+Ky+PoUBRcgUhFLKLiJ3oHknmYHfK6V2i8hK13bfeYdu24ZK1r4kOTnZ78PGjb95lGC9TKxWK0uXLgVg5syZvPPOO373O/fcc/nwww/ZsGED3/ve9/jXv/6FUsqvAps6dSo33HADy5cvZ/ny5X6P9+6773qZsOrq6qivr2fDhg28/LJWyvrSSy8lPT29Q9tPP/2U8847D7eZ8JprruHAgY4KPRA5AK6++mrMZrPfbV/5ylcwmUwUFBQwevRo9u3b1+lxumLjxo0su+JqTCYz2YMGM2/+Ag7s2kFWRhpz5swhJycHgMLCQoqLiykqKgK0+5kcF4PFZKKuRbOpOxXUNNlIT7RiiiCvot7E5rTx1tG3qGmtAcAkJi4cdSH5afldtjOJiUlZk0iyJvGvo//CoRw02Zt47fBrfGXcV0iMSewL8QckIQ2UU0q9Cbzps86vYlBK3dxd22iguxHEhAkTqKmpwW63Y7FYKCkpYdiwYR32tVgsOJ3tZgm9H3tMTIxHqZjNZux2e4f2AAsWLGDjxo0cO3aMyy+/nMceewwR8SgXPW+88QYbNmxg3bp1/OhHP2L37o762ul08vHHHxMf39HbJBAlF8g+gcgBkJjY+UPDtx8R6fJ6+kMpRYvNgVNnLbOaTVgt2gggNjbWs76zexBvNWM2QXWTpiTsTkVds43U+JiIcj3tDZzKyX+O/YdTjac86y4YeUG3ykHPqJRRXJZ/GW8cfYM2RxvN9mbeOfYOl+Vfhkmie+QVLqI+klpPZ2agUJKcnEx9fb3XclcjCIDFixezdu1arr32Wp5//nkuv/zyDvvk5uby1FNP4XQ6OXnyJJs3bw5aloULF/L973+fhQsXYjKZyMjI4M033+TRRx/1aud0Ojlx4gSLFy+mqKiIF154gYaGBpKTk6mrq/Pst2TJEp588knPXMv27dspLCxk4cKF/OUvf+H73/8+b731FtXV1R1kmzt3Lv/7v/9LVVUVKSkp/P3vf2fatGk9kqM7/v73v7NixQqOHj3KkSNHGDduHPX19Z1ez5iYGGw2m1cqioZWO3PmzedPf/gd11z/VRzN9WzatJGf/Wx1UCMSq8VMcpyivkVTIK12J41tDpJio+unuaNiB0dq2+cczh12LmPSxwR9nKFJQ7k492JeO/waCkVpQymbyzZzztBzelNcAxeG2g0xmZmZzJ8/n8mTJwc0SQ3w2GOP8fOf/5wxY8ZQVVXF17/+9Q77zJ8/n7y8PKZMmcJ3vvMdZsyY0e1xr732Wn76058yffp0Dh8+jLu40sKFCwEoKioiLS2tgwnI4XDw1a9+lSlTpjB9+nTuuusu0tLSWLZsGa+88opncviJJ55g69atTJ06lYkTJ/L009pg8cEHH2TDhg3MmDGDt99+m5EjR3aQbejQoTz00EPMmzePCy64wO/5BCpHd4wbN45FixZxySWX8PTTTxMXF9fl9bzttts8pi03TW0OvrTsciZMmswFRXP40kUX8pOf/IQhQ4Z0278v7nkJN42tdtrs0TNpfablDJtPtSvcqdlTKRxU2OPj5STnMHtIe9jUZ6c/41hdpyUNDM6CqKpJPWvWLOVbMGjv3r1MmDAhTBIZRBtOpTjT0OaJcYi1mHrFJKSU5iZrc3kzWUxCRqK1z0xNofqdOJWTVw6+wumm0wBkxWdxZcGVmE3+54cCRSnF60de50S9Fi4VZ4njuvHXGcF0PUBEtimlZvnbZowgDAyCoFEXACcCyXG9M18gIqTGWzyxEnan6hBT0R/ZUbHDoxxMYuL8keeftXIA7Xr918j/8kxQt9hbvEYpBr2DoSAMDALE5nB6PbSTYy2YezFdhtlkIlE399DQau/X8RHVLdVeD+1Zg2f1qltqQkwCC3MWepb3VO2hoinysin0ZwwFYWAQAO40Gm5izCbiYs7+TdiXBKvZK5CuvsXeb9NxfFz6MQ5XkoSs+CymD+oY13K25KbkMjJFm9NSKDaUbOi31ysSMRSEgUEAtNic2BztD56UOEtI5gdEhJT4dm+pNoeT1n44YV3aUEpxXTGgpc9YPGJxr5iWfBERioYVedxcTzedZn/1/l7vZ6BiKAgDg25QStHY2h7HkGg1YzGH7qcTYzYRb/X2aupPb8VKKT459YlnuSC9gOyE0OVJS4tL8/KK+rj0Y9ocvZPgcqBjKAgDg25otjm8JqYT+iBGIdHqPWHd0o9GEe603aBNTM8ZOifkfc4cNJOkGC1HVrO9mZ2VO0Pe50DAUBAhJlTpvrdv386bbwYeaL5161b+53/+p9PtxcXFvPDCC0HJCYGlzy4uLmby5Mm9Ikdpaaknl5Rvgr1AWLNmDaWl7Xkfb7311i4z3Gqjh/aJ6USrpU9SYZhNQkI/HEU4ldNr9DAlawop1pSQ9xtjjvGKjdhRscNIDd4LGAoixIQq3XewCmLWrFk88cQTnW7vqYIIlrORw263M2zYsKDrOejxVRDPPfccEyd2Xqywqc3hScdtErxMP6EmwWrBrYscTi21R6Sz78w+T64lq9nKjMHdB3D2FmPTx3qUUYu9xRhF9AKGgggxoUj33dbWxgMPPMBLL71EYWEhL730ElOmTKGmpgalFJmZmZ4COzfeeCPvvvsu69ev9+RY+uCDDygsLKSwsJDp06dTX1/Pfffdx8aNGyksLOTxxx/H4XBw9913M3v2bKZOncozzzzjke+OO+5g4sSJXHrppZSXl/s9j23btjFt2jTmzZvHr3/9a8/6YOVYs2YNV199NcuWLWPJkiVeoxGAEydOcPHFFzNu3DgefvhhgA77rF69moceeoi1a9eydetWbrjhBgoLC2lubua8887DHVz517/+lSlTpjB58mTuvfdenErR1GZn9LAsHv3hg5w/fw7nzpvH6dOnu72PvYHJJCRa9W6vjogeRTiVk8/LP/csF2YX9mngmtlk9lJI2yu2G6OIsyS6Er50x/uPdr9PT1n8Xb+rQ5Hu22q18sMf/pCtW7fy5JNPAvD+++/z4YcfMmrUKEaPHs3GjRu56aab+OSTT/jNb36DPsJ89erV/PrXv2b+/Pk0NDQQFxfHqlWrWL16tac+wrPPPktqaipbtmyhtbWV+fPns2TJEj7//HP279/Pzp07OX36NBMnTuSWW27pcC5f+9rX+NWvfsWiRYs6VYyByLFmzRo+/vhjvvjiCzIyMiguLvY6xubNm9m1axcJCQnMnj2bSy+9lKws/772V111FU8++SSrV69m1izvwNHS0lLuvfdetm3bRnp6OkuWLOFva19m8UWX0tTYyKw5c3n8p6u49957PenT+4J4q5nGNjtKaVHczTYHCdbI/NkW1xZ76jRYzVamZk/tcxnGpY9j2+lt1LfV02JvYVfVrpC41w4UjBFEH+NO1ufvb+LEiX7fEANxp1ywYAEbNmxgw4YN3H777ezcuZOTJ0+SkZHhKXDjZv78+Xz729/miSeeoKamxqOM9Lz99tv88Y9/pLCwkLlz51JVVcXBgwfZsGED1113HWazmWHDhnH++ed3aFtbW0tNTQ2LFi0CtFGMPwKRA+DCCy8kIyOj022ZmZnEx8dzxRVXsGnTpi6vU2ds2bLFk27cYrFw/fXXs/6DDYCmkL98+TJEhJkzZ3ZQUqHEJN6jiKa2yBxFKKX4rPwzz/LkrMlYzdY+l8N3FPF5+efGKOIsMBREH1NfX+8xq/j+7dmzh6ysLE+6b6DTdN++LFy4kI0bN7Jx40bPg27t2rV+Ryv33Xcfzz33HM3NzZxzzjl+s48qpfjVr37lUV5Hjx5lyZIlQPcKSykVkFILRA4Ifeput8x67E6Fu8qtJSaGeNdDuqv06aEiPsbs8WhyOFVEJvIrbSylvEkzN5rFzNSsvh89uBmfPt7j0dRib+FAdegKhUU7IR2risjFwC/Riv48p5Ra5bP9cuBHgBOwA/9HKbXJta0YqAccgL2zZFJB0YkZKJSEKt2373FHjBhBZWUlbW1tjB49mqKiIlavXu0xQek5fPgwU6ZMYcqUKXz88cfs27ePESNGeB3voosu4je/+Q3nn38+MTExHDhwgOHDh7Nw4UKeeeYZbrrpJsrLy3n//fe5/nrvYn9paWmkpqayadMmioqK+Mtf/uL3PAORozveeecdzpw5Q3x8PK+++iq///3vGTx4MOXl5VRVVZGUlMTrr7/OxRdf7Pe6uXGnG6+srCQtLY2//vWv3HLb7YBWVzqcRXxMJiHeavak+WhqcxAbgijus0E/9zAuYxwJMQld7B1azCYz0wZN48OTWtnaLyq/YGLmxKirsdEXhGwEISJm4NfAJcBE4DoR8XUX+Q8wTSlVCNwCPOezfbFSqrBXlEOYCFW678WLF7Nnzx7PJDVoD7mxY7WaFwsWLODkyZOeKmZ6fvGLXzB58mSmTZtGfHw8l1xyCVOnTsVisTBt2jQef/xxbr31ViZOnMiMGTOYPHky3/jGN7Db7Xz5y1+moKCAKVOmcPvtt3vMSL784Q9/4Fvf+hbz5s3zW0AoUDm6o6ioiBtvvJHCwkKuvPJKZs2aRUxMDA888ABz585l6dKljB8/3rP/zTffzMqVKz2T1G6GDh3Ko48+yuLFi5lWWMjkqYVcfOmybvvvK/TeU20OpyfrayRQ2VzJ8TqtnrkgFGb3PJV3bzEhY4LHxFXdUk1JfUmYJeqfhCzdt4jMAx5SSl3kWv4ugFLK70yxa//fK6UmuJaLgVlKqcpA+zTSfRv0BjVNbZ70FvExZq/UF+FEL1dcjJnUXpTrbH4n7x1/j31nNPPg6NTRXJx3ca/JdTZsLNnocXUdmTKSpaM7Vko0CF+67+HACd1yiWudFyLyZRHZB7yBNopwo4C3RWSbiNzWWScicpuIbBWRrf7qORsYBIPdJ/dRX8Y9dIfee6nF5sDhDP9kdYu9hUM1hzzLZ1MIqLeZmj0Vcc3eHK87Tk1LTZgl6n+EUkH4M/h1+EYrpV5RSo0HlqPNR7iZr5SagWai+paILPRt62r/rFJqllJqlrvgvYFBT2nWBaNZLSZiQphzKVh85Wlu69vJcn/sr96P3anJkRWfxeCEwWGWqJ3U2FRPplfQ5iIMgiOU3/4SYIRuOQco7WRflFIbgHwRyXItl7r+lwOvAKFP6GIwoFGuOAM3CRE2EQx4pd9otjnD6vKqlGJ35W7P8qSsSRE3Eaz3ptp/Zj+tjtYwStP/CKWC2AIUiEieiFiBa4F1+h1EZIy4vlEiMgOwAlUikigiya71icASYFcIZTUwoMXmxP28NYlgtUTO6MFNrMWEu1yEU6mwpgI/2XDSK63G2LSxYZOlM3KSc0iP02qs25w2DlYfDLNE/YuQ/QKUUnbgDuDfwF7gb0qp3SKyUkRWuna7EtglItvRPJ6uUdor0WBgk4jsADYDbyil/hUqWQ0MwNu8FG81R9zbMGhxHvpCReHMz7Srqv2dbWz6WGLMkTGZr0dEmJzZnnZlT1XniRkNOhLSOAil1JvAmz7rntZ9fgx4zE+7I8C0UMpmYKDH7uM6Gh+B5iU38THtMRGtdicOpxOzqW9HO422Ro7WHvUsT8qc1Kf9B0NBegEflX6EQzmobK6koqkipPUpoonIG0NHGT3J5vree+954g9WrFjRaeTu559/zq233gr0LPU1dMyeunPnTm6++eZO97/uuuuYOnUqjz/+OA888ADvvvtu0H2Cd9I+PcFmqT0b9Jld9aOHWIvJq9Z0pKUrLz9dhtU1Wf3tO27n8x19b33dW7XXM/8xNHEomfGZfS5DoMRZ4shPy/cs7z2zN4zS9C8MBRFiglUQTqeTFStW8OKLL7Jr1y5GjRrF888/73ffH//4x9x5551nJZ/vA23KlCmUlJRw/PjxDvuWlZXx0Ucf8cUXX3DXXXfxwx/+kAsuuOCs+velJwqip6kv3ArCd3K6q9FDpKQrd7vf/vzJ3zBqzLg+naxWSnniHkCbnI50JmS2x3gcqD6AzWnkZwoEQ0GEmGDTfVdVVREbG3ZA0WcAACAASURBVOuJiL7wwgv5xz/+0WG/+vp6vvjiC6ZN62iJe+2115g7dy7Tp0/nggsu8KSnDiS9NsCyZct48cUXOxx3yZIllJeXU1hYyMaNG72KBeXm5vLggw8yY8YMpkyZ4smrtHnzZs4991ymT5/Oueeey/79ndcL9pfGvLP2vmnAnU4n3/zmN5k0aRJLly7lS1/6kke2bdu2sWjRImbOnMlFF13EqVOnvFJ/TyucTlOTFlXtnpyO9HTl58yeSWtLM1++dAmfbdtKq93ZIV25m6SkJO6//36mTZvGOeecc9bpyksbS6lrqwMg1hzL6NTRZ3W8vmBY4jBSY1MBaHO0caTmSJgl6h9EZt7gEPHU9uBMPcHwzcJv+l0fbLrvCRMmYLPZ2Lp1K7NmzWLt2rWcOHGiw75bt271eojoKSoq4pNPPkFEeO655/jJT37Cz372s4DSa4P2lrxq1Sruuecer+OuW7eOpUuXes7Ft5BRVlYWn332GU899RSrV6/mueeeY/z48WzYsAGLxcK7777L9773Pb8KD/ynMa+rq+u0vT4N+Nq1aykuLmbnzp2Ul5czYcIEbrnlFmw2G3feeSf//Oc/yc7O5qWXXuL+++/n97//vSf1d/6EqbS55h/ck9P9IV15fUv7W/DR4yc6pCt/9dVXWb58OY2NjZxzzjk88sgj3HPPPWedrlw/eihIL8BiivzHiIgwIWOCp9rdnqo9jMsYF2apIp/Iv7NRRiDJ+l588UXuuusuWltbWbJkid802KdOnaKzwMCSkhKuueYaTp06RVtbG3l5eUB7eu0bbriBK664gpycHL/tBw0a5FV1LVCuuOIKAGbOnMnLL78MaKm/V6xYwcGDBxERbLbghvZdtdenAd+0aRNXX301JpOJIUOGsHjxYgD279/Prl27uPDCCwFwOBwMHTrUcwyH0+lRDgDxMSa/6crfeuutDrIFej0DSVcOeNKVL1++PKBrozeFbd68hUWLzvN8J2644QY2bNjA8uXLsVqtnlHPzJkzeeeddwI6vj9sDpvX2/e49P7zkB2fMZ5Pyz5FKcWpxlNUt1R7XGAN/GOYmPqY7tJ9A8ybN4+NGzeyefNmFi5cSEFBQYfjxMfHd5rG+s477+SOO+5g586dPPPMM579Ak2v3dLS0mmCva6IjY0FvFNi/+AHP2Dx4sXs2rWL1157LaDU23q6aq9PA96ZDV4pxaRJkzxpy3fu3Mnbb7/t2a5PnW01mzCbTP0mXbnFbPKkK1BK4VD+YyJiYmI8/ZxtuvIjtUc89vv0uHQGJQzq8bH6moSYBHJTcj3L+6s7N3caaAyoEURnZqBQ0pN03+Xl5QwaNIjW1lYee+wx7r///g77TJgwgZ/97Gd+29fW1jJ8uJb2Sj/BHWh67QMHDnRqvgoWvSxr1qzpdn/f6xVo+6KiIp5//nlWrFhBRUUF69ev5/rrr2fcuHFUVFTw8ccfM2/ePGw2GwcOHGDSpEkkJydTVVPrOYY7vqA/pSt3e1vNmDWbH9z3HSorK0lPT+evf/3rWTsw+ENvXhqfMT4iY0W6YlzGOI977sHqg8wdMrffnUNfYowgQkxP0n3/9Kc/ZcKECUydOpVly5b5rdo2fvx4amtr/T6MHnroIa6++moWLFjgZc8ONL32+++/z6WXXtrDM/bmnnvu4bvf/S7z58/H4eg+qMs3jXmg7a+88kpycnI8qcnnzp1LamoqVquVtWvXcu+99zJt2jQKCwv56KOPAPjqTTfxnf+9k/8qmktLczOxMe0/h/6Srtxdp2LwkKF878Efct7ixUybNo0ZM2b4rSNyNtS11XGyQSt/Kwhj0yMvcro7RiWPItasjXTr2+o51XgqzBJFNiFL9x0OBlq678cff5zk5GRPLERv0NrayqJFi9i0aVOnJUAjlYaGBpKSkqiqqmLOnDl8+OGHDBkypNP961tsnoCz3k6f3Zfo04AnWs0kxQV/HoH8TraUbWFL2RYARiSPYFl+5NTLCIYPTnzA7ioth9TEzImcN+K88AoUZsKV7tsgxNx+++0eu39vcfz4cVatWtXvlAPA0qVLKSwsZMGCBfzgBz/oUjkopbzSVMTF9N+fgj71RqgS+CmlvEp3js8Y38XekY1+5HOo5pAnG61BR/rfU8DAQ1xcHDfeeGOvHrOgoMDvpHh/YP369QHv22Z34tQn5ougtN7BEmsxIQJKaQn8bA4nVkvvpgqpaK6gtlWbr7GareSm5vbq8fuSIYlDSLGmUNdWR5ujjeN1xxmdFvmxHOGg//4qgiCazGgGvUOzz+ihP09UighxFn0Cv+AyvAby+9BnQc1LzSPG1D/NcaBdL/0owvBm6pyoVxBxcXFUVVUZSsLAg1MpL/fWuAhOzBcoXhle7Y6Av+9KKaqqqoiLi+t0H6dyelWNK0jrnyNMPXoFcazuGC324NyvBwpRb2LKycmhpKQEoxypgRubw0mza3LaZBLOxEbHz6ChxY7TpRhOW81YAjSbxcXFdRrkB1DaUEqjrRGAeEs8Ocmd79tfSItLY1DCIMqbynEqJ4drDveLnFJ9TXT8MrogJibGE0lsYADwz+0nOVKlPfDOzc9kwujIzUQaDJsOVrKl+AwAE4Ymc/GEod20CAz96CE/LR+TRIfhYWz6WMqbygE4WHPQUBB+iI47bWAQIC02B8eqmjzL44Ykh1Ga3mXskCTP58MVjV71LXqK3WnncM3h9j76YexDZ+Sn5SOuWPRTDac8oySDdkKqIETkYhHZLyKHROQ+P9svF5EvRGS7iGwVkaJA2xoY9IRD5Q04XO5Lg1PiSEuwhlmi3iM7KZaMRO182uxOiivP/oF3ov6Ep45zsjWZwQmDz/qYkUJiTCJDk7RRlkJ5KUIDjZApCBExo5URvQSYCFwnIhN9dvsPME0pVQjcAjwXRFsDg6DZX9YeeR5NowdweecMbj+nfWWBp/zoDL330pi0Mf3a28sf+gl3Q0F0JJQjiDnAIaXUEaVUG/Ai4BX7r5RqUO3uFomACrStgUGwNLbaOVHdbl4aOzipi737J3qlV1zZSKu95zWrbU4bxXXFnuWC9P7vveRLXmqeR+mdajxFfdvZK9VoIpQKYjigL2RQ4lrnhYh8WUT2AW+gjSICbutqf5vLPLXV8FQy6IqD5Q24X0eGp8eT3IOUFJFORqKV7GQtut7uVBwu77mZ6XjdcU+UcVpsGplx0TGZrychJoGcpHavLGMU4U0oFYS/sWgH52yl1CtKqfHAcuBHwbR1tX9WKTVLKTWrs/oIBgYAB07rzEuDo8u8pEc/itCfc7DovZei0bzkRl+vWn/OBqFVECXACN1yDtBpFRql1AYgX0Sygm1rYNAdDa12Smu0sqIiMGZQ9JmX3OjnIY6fafLKORUoNqeNY3XHPMv6h2i0MTp1tMd1t7yp3JNSxCC0CmILUCAieSJiBa4F1ul3EJEx4notEZEZgBWoCqStgUEwHNKZl3LSE0iMkuA4f6TGxzA4RYuMdjgVRyqCNzP5mpcy4vxXxIsG4ixxjEhufx81zEzthExBKKXswB3Av4G9wN+UUrtFZKWIrHTtdiWwS0S2o3ktXaM0/LYNlawG0c9BnamlIIpHD24KdBPwB8uDNzMNFPOSmzFpYzyfD9caCsJNSF+jlFJvAm/6rHta9/kx4LFA2xoY9ITGVjsnB4h5yU3BoCQ2HawE4FiVZmYKNOfUQDIvuclNzcUkJpzKSUWTlrk2NTY13GKFHSOS2iDq0ZuXhqfFR7V5yU1agpVBKZo3k8OpOBpE0NxAMi+5iTXHeuWYcpclHegYCsIg6jlY3uD5XBDF3ku+FAxqP1f9NeiOgWZecpOf2j5SMuYhNAwFYRDVNLXZKXEFxw0U85Ib/VzLsQCD5mxOG8frjnuWB4J5yY0+aO5002kjaA5DQRhEOXrz0rDUeJIGgHnJTbpP0FwgZqYT9SewOW3AwDEvuYmzxHkFzR2pPRJGaSIDQ0EYRDUHT+vNSwNn9OBGP4rQX4vOOFLT/lAcnTZ6wJiX3OhLj+qvxUDFUBAGUUuLzUFJdbNneSCZl9zo51yOVTV6VdLzxeF0eOVeGp068Oo056XkeVKAlzWWDfgU4IaCMIhaDlc0eCqsDU2Ni8rcS92RkWglM0lLAW5zKI5Vdf7AO9lwkjZHG6Cl9s6OH3ipaxJiEhiWNAzQUoAP9FGEoSAMopZDOs+dgTh6cDMmu/3cD3XhzaT33BmdOvDMS270E/MDfR7CUBAGUUmr3cFxXeW4Aa0gdHMvRyobsfupNOdUTo7Wtfv+DyTvJV/yUtvNTKWNpTTbm7tpEb0YCsIgKimubMLuqhyXnRwbVZXjgiU7KZbUeM281mZ3cqK64wPvVOMpWuwtgFZpLZoqxwVLYkwigxO181dKUVxbHF6BwoihIAyiEsO81I6IeF0Df2YmvXlJHw8wUNFP0A9kM5OhIAyiDpvDSbFuMnagKwjwvgaHKxpwOtvLqyilvFJLDETvJV/yUvM8n0/Un/BM3g80DAVhEHUcP9PkcedMT4ghM3HgmpfcDE2N8wQJNrc5PMkLQYsadrtzxlniPF48A5nU2FSy4rMAbX5Gn7xwIDFwwkoNBgyHfHIvDUhzia0ZGiuhtQ5aahFbE3McNRyprkdhomLPIEYUjITYFI7Ut5tQclNyPcVzBjp5qXlUNmsZcY/WHo3KmtzdYSgIg6jCt0BOfvYAMS/Z2+DMYag+BrUnNOXgQ26zjab6Om3hsAllSwMFRxsOgNkCscmMTp0IDru2PMAZnTqaLWVbADhWdwy7047FNLCuy8A6W4Oo52R1s6fEZnKchcGulNdRidMJZ47A6V1QdVB7sHdBSpwFi0mwOxWtdieNrQ5aLXZqHa3gaCXG1kLO0Q/hxFbIKoAhUyA9T8tyOADJiMsgNTaV2tZabE4bJfUl5KbmhlusPiWkCkJELgZ+CZiB55RSq3y23wDc61psAG5XSu1wbSsG6gEHYFdKzQqlrAbRwaGK9gyc+YOSotO85LBB2U4o2QJNZ/zvIyZIyoa4NIhLAWsSIiYcljMcr6gnxtFIukWotxzVFIBSjIxJxiIm7fin92h/iVmQMwsGTwbzwIpEFxFGp47m8/LPAa3SnKEgegkRMaOVEb0QKAG2iMg6pdQe3W5HgUVKqWoRuQR4Fpir275YKdVxrGxg4AelFIfLdd5L0WZecjqgdDsc2wRtTR23J2ZB9jhIGwnJw8DScXI+O66Bj1tLAaiLjyE+Mw5SMqC1gdFJudBYB83V7Q0aK2H/v6B4E+QugCFTwTRw5ij0CuJY3TGcyjmg5mhCOYKYAxxSSh0BEJEXgcsBj4JQSn2k2/8TIAcDgx5SVtdCQ6tmZom3mhmeFh9miXoJpaDqEBx+H5qqvLdZYmFYofaGn5jdrTloZEYCVouJNruTsoZqLDGnSbBaMCVkMGriVWCKgYZyKPsCTu3QRhMArQ2w/y1t1JJ/PmQOjEjrQQmDSIxJpNHWSIu9hdKGUq/Kc9FOKBXEcOCEbrkE79GBL18H3tItK+BtEVHAM0qpZ/01EpHbgNsARo4ceVYCG/Rv9N5Lo7MSMZmiwLzUWg8H/g2VB73Xx6VAzhwYOlVTEgESYzYxKjOBg6cbqLWfRBrbSLBayEnOwWp2jTiSB0PyhdqI4dQOOPEptLlGZo2V8MXfYPBEGHMBWBN76UQjExEhLzWPXZW7AC1obiApiFCOlfz9OpWfdYjIYjQFca9u9Xyl1AzgEuBbIrLQX1ul1LNKqVlKqVnZ2QMv+6SBhlIquqKnlYJTX8Dm33orB4sV8hfDnG/AiNlBKQc37mtTYy/hTKMWAOY3OC4mDkbOhbkrIW+B9xzE6T2abKd3B91/f0N/bY7WHkUpv4+xqCSUI4gSYIRuOQco9d1JRKYCzwGXKKU842elVKnrf7mIvIJmstoQQnkN+jFVjW3UNGnmEKvFxMiMhDBLdBbYWmD/G1BxwHv9sOmQtxCsZ3duuZmJOGmh0VGFsmseTbkpuZ03sFght0jr//B7ULbLJWcz7FmneVIVXOR3ziMaGJY0jFhzLK2OVhptjZQ3lXtyNUU7oRxBbAEKRCRPRKzAtcA6/Q4iMhJ4GbhRKXVAtz5RRJLdn4ElwK4QymrQz9GPHnIzE7GY++lEYt0p2PYHb+UQnwaF18G4i89aOQDExZiJS6jyvAmLI52EmACOa02ECctg6lc0E5ebsl2wbQ3Unz5r2SIRk5i8vJcGUm6mkP2KlFJ24A7g38Be4G9Kqd0islJEVrp2ewDIBJ4Ske0istW1fjCwSUR2AJuBN5RS/wqVrAb9n8MVUWBeKv0cPv8TNNe0rxs+A2Z9HdJze7Urk7XC89nZOii4xpn5MPtWLU7CTVMVfPbHqDU5+ZqZBgohjYNQSr0JvOmz7mnd51uBW/20OwJMC6VsBtFDbZON8rpWAMwmITern5mXnE448h6c2NK+zmKFcV+CQRN6vbtWRyt2U6U7/AFbSxZNbXYSrEE8DiyxMGEppI/SJtEdNnDaNZNTY6VmCouiGJSc5BxiTDHYnDZqWms403KGjLiMcIsVcvrpONzAoJ1DutHDyIwEYi3mMEoTJPZW2PUPb+WQNAhmfi0kygE0f36zGZJjLSSY07FKold6kqAYMkWTNSFT18FHsPvldhfZKCDGFMOIlPYp1YFSitRQEAb9nsP91XuprRG2/0WLcXCTVQDTb4SE0L2dum3o6YlW0izDga5LkXZLYibMuAkydJ5QFQdgx1+1iewoYSDWiDAUhEG/prHVTmmt9hASgdHZ/cQvv7kGPv+z98TuyHNg8pUh9QayOW0crzsOQEailVSL5tN//EwTrXZHzw8cEwdTrtZcb93UntTOsaXubESOGEaljPJEUVc2V1LXFh3n1RUBKQgRWSoygOLLDfoNRyoacbulD0uLD86OHi4aK7UHpzuPkgiMu0SLbwix3b6kvgS7U4s2H5KUycg0LXbI4VQUV/pJ3xEMJpMWPDfmgvZ1jZXaxHtnOaP6EbHmWIYnDfcsD4TJ6kAf+tcCB0XkJyISGsOogUEP0Cfn6xfmpYYKTTm0uuQ2mWHSl7V0GX2A3jSSl5rHmEHJnuWzMjPpGTFbc4d1v1O21GmmtChQEqPTdGamATAPEZCCUEp9FZgOHAb+ICIfi8ht7lgFA4Nw0GJzcOJMu4074ms/NJRrD0q3Xd4co8UUZI/rk+6dyklxbbFneXTqaC+lWlzViN3h7J3OhkyGKVe115VobdDOvbGq63YRTl5KHuJKElHWWEaT7SxHXRFOwGYjpVQd8A/gRWAo8GXgMxG5M0SyGRh0SXFVIw5XbeXBKXGkxkdwOuqGctj+QrtysFhh2nW9Ht/QFaUNpbQ6NHfgxJhEBiUMIjPRSnqCdt3a7E6On+nFB15mvjYvEUVKIiEmgSGJQwBQKIrrisMrUIgJdA7iMle6i/eAGGCOUuoStFiF74RQPgODTtGbRPIjeXK66Yy3R4/FClOvhdThXbfrZfTmpdGpoxERRIR83Sii18xMbtJzYcpX2pVEW6N2LfTBgP2MvNQ8z+do92YKdARxFfC4UmqqUuqnSqlyAKVUE3BLyKQzMOgEm8NJcaWu9kOkzj+01GoPRHf9hjApB6WU16Sq3pauv3ZHKhtxOns5GV36KJh6jW4kUQ87Xmyfh+ln6K9dSX2JZ1QWjQSqIE4ppbwS5YnIYwBKqf/0ulQGBt1wrKoJm0N7kGUkWslMisDSoq0N2oPQ7eZptmgPyj5WDgCnm07TaNMUarwlnqGJQz3bhqTEkRynPbyb2xycrAlB7ELaSJh8lTYpD1pRoh0v+i98FOGkWFPIis8CtHmdY3XHwixR6AhUQVzoZ90lvSmIgUEwRHxqb3srfPFSu+eOyazFOKSGp5aA3uMmNyXXqyqaiHhN8B8sD9GbfUYeTFze7t3UWAm71vbLiOuBEjTXpYIQkdtFZCcwXkS+0P0dBb7oGxENDLxxOBVHKtsVREGkKQinA3a9rE1MgxbbMPFy70jjPkQp5T3/kNZRDr2SPVzeGLqaB9ljtRxO7niP2pOw559aPqp+hP4anqg7gc3Z/5RcIHQ3gngBWAb80/Xf/TfT5fpqYNDnnDjTRKtNe6CkxMeQnRxB5iWlYN8bUF3cvm7cJX3myuqPqpYqT9Sv1Wz1CvZyMzwtnnirZv5paLVzqrYldAINnuQdTFd5EA7+G/pRIZ702HTSYtMALTr9RP2Jblr0T7pTEEopVQx8C6jX/SEi0Z/K0CAi8TUvSSRlDT2y3jvldd4CGBrexMR689KolFFYTB2jzU0mXzNTL3sz+ZIzS6tW56Z0u5bkr58gIgMiaC6QEQTANmCr6/823bKBQZ/idCqv2g8RZV4q3Q7HP2lfHlYIo+aHTx4Xh2sPez77LS3qosDH3TXkpTVHL9ZGE26ObtBKmfYT8lPzPZ+L64pxOM8il1WE0qWCUEotdf3PU0qNdv13/4XHoGowoDlZ00xTm/ZDTIq1MDQ1LswSuagu1uoiuMkco5XhDPPoprqlmuqWagAsJgsjk0d2uu+IjARiY7RHQl2zjYr6ELtvisD4SzU3WDf73oDaktD220tkxWeRbNWSSbQ52ihp6B9yB0OggXLzXaU/EZGvisjPXeVCu2t3sYjsF5FDInKfn+036Ca+PxKRaYG2NRiY6Gs/5A9KjAzzUmOVNimtXBOtSYO0SWlT+PNbHq5pHz2MSB5BjLnzaHOzSRid1R5wGHIzE7TnonLXk3DatfoYzdWh7/ssERGvUUQ0mpkC/Qb/BmhyPcDvAY4Bf+qqgYiYgV+jucNOBK4TkYk+ux0FFimlpgI/Ap4Noq3BAEMp5VX7oWBQBKQCszXDzr9rbq0AsUlaeokQpuwOBr33Un5afhd7aoQkeV93xMTD1Ku1/6DFRuxc235NIxj9PMTRuqM4Vf/yxuqOQBWEXWkGycuBXyqlfgl09+ucAxxSSh1RSrWh5XC6XL+DUuojpZT7VeETICfQtgYDj7K6FupbtFTV8VYzw9PiwyuQ0wm7X21/2zVbtGCwuJTwyuWitrWWyuZKAExiIjclt9s2ozITiDFro7IzjW1UNvTRQzo+XUvu555Ab6yEva9FvGfT4ITBJMZoo64WewulDaVhlqh3CVRB1IvId4GvAm+43vC7y4w2HND7fpW41nXG14G3gm3ryiq7VUS2VlRU+NvFIEo4cFqfeykJkynM5qXD73m7s45fBilDO929r9Gbl0Ymj8Rq7n5UE2M2kZel82Y63UejCNCCCMdd3L5ceRCOftB3/fcAEYnqoLlAFcQ1QCvwdaVUGdrD+qfdtPH36/X7OiAii9EUxL3BtlVKPauUmqWUmpWdnd2NSAb9FaUUB0+3R/iG3Xvp1A4o0dWRzi2CQePDJ48f9N5LgZiX3BQM7oOo6s4YMgVGzGlfPvaxt9twBOLr7hpNZqZA60GUKaV+rpTa6Fo+rpT6YzfNSoARuuUcoMP4S0SmAs8BlyulqoJpazBwOF3X6jEvxcWYGZGRED5hak96eyxlj9MURARR11ZHRZM2ojaJidzU3IDb5mYmesxMVQ1tVPWVmcnN6MVaqnA3+9/0Ls0aYQxNHEqCRfs+NtmbKGssC7NEvUegXkxXiMhBEakVkToRqReR7gqybgEKRCRPRKxoVenW+Rx3JPAycKNS6kAwbQ0GFgd0o4f87ETM4TIvtTbA7pe1dBoASdkwfmnY3Vl90XvU5CTnEGsOPNrcajGR29feTHpMJphwWbtnk8Pl2RShif1MYvIaRRyqORRGaXqXQE1MPwEuU0qlKqVSlFLJSqkuZ+KUUnbgDuDfwF7gb0qp3SKyUkRWunZ7AMgEnhKR7SKytau2QZ+dQVSglPJ6SBUMDpP3ktMBu1/RlARATJyWgC9CPJb06Ocf9K6YgaL3ENOb9voM32vbUhvROZu85iGiyMwUaIX300qpvcEeXCn1JvCmz7qndZ9vBW4NtK3BwOR0XSt1zVoytNgYEyPDZV469J/2IC53Ar749PDI0gX1bfWcbtJMMiLiVeAmUPKyErGYBLtTUdnQxpnGNjIS+1gRJmZqE/+7/qEtVxfD0fWQf37fyhEAw5KGEW+Jp9neTJO9iVONp/zmvOpvBDqC2CoiL4nIdS5z0xUickVIJTMwcKGfKM3PTgqPealsJ5zc1r6ctyhs2Vm7Qz96yEnKIc4SfLS5r5npQDhGEaBlf83VpSs5/imU7wuPLF1gEpOXI4D+HvRnAlUQKUATsIT2jK5LQyWUgYEbzXspzLmXGsrhwL/al7PHwchz+l6OANHbwIPxXvLF25upj+ch9OQu0FKXuNn/RkTWtdZf62gxMwVkYlJKfS3UghgY+KO8vpVanXlpVGYf1562tWhpNByaBxUJmVr+oAiblHZT21pLeZNWh8LXRz9YvMxM9a3hMTOBdq0nLIVta7Ra1vY2zVFgxoqImv9xezM12ZuixswUqBfTWBH5j4jsci1PFZHvh1Y0AwPYXxZG85JSsO91XaR0DEy+AiwRVH/CB9/cSz0xL7mJtZi9zEz6e9HnxMTDpCu8I633vxlRkda+3kzRYGYK1MT0W+C7gA1AKfUFmuupgUHIUEp52b7H9bX30vFPtGheN+OXQmJW38oQJHrzUkFawVkfb9yQ9mt+4HR96FOAd0XyYBh7Ufty+V44+Vn45PFDtJmZAlUQCUqpzT7r7L0tjIGBntJa79xLfRocV33MO83DiNkRFyntS01LjXfupSCC4zojNzMRq0V7TGi5mdrO+phnxdCpWp0NN4f/owUuRgi+QXMnGyJHtp4QqIKoFJF8XOkuROQq4FTIpDIwAA6UeafW6DPzUmu95nPvfltOzdGieyMc/ehhVMqooILjOsNqMXmlAA+bN5OeMRdqownQYlP2vBoxQXS+3kyHqvt30FygCuJbwDPAeBE5CfwfYGXXTQwMQgOIcQAAIABJREFUeo7T6W1eGttX5iWnU1MObY3asjUBJi3X6hZEOL3lveTLWJ2ZaX9ZmM1MoGXNnfTl9rmglrqIyvxakN5u2jtcexi7s/8aW7pUECLybRH5NrAcLWjtEeBptPQYV4ZePIOBSkl1e+W4xNg+TO199AOocSUSFtFSPsRGQN2JbjjTcoYzLWcArXJcXkrwwXGdMUpXaa622cbpugio0xCfDhOWtS+fORIxNa0HJwz2qjR3ov5ENy0il+5GEMmuv1nA7UA6kIY2ejAK+BiEjP36zK2Dk/smtXflIe+a0rkLIKP3HrSh5GB1+2T6yJSRXVaOCxaL2cSY7PaYiP2RYGYCyCrwjkcp3uidfj1MiAhj0trjNvT3pr/RXU3qh5VSDwNZwAyl1HeUUv8XmEl7cR8Dg17F4VRe1cz6xHupuQb2vda+nDEaRp0b+n57AaWU10NobNrYXu9D7810MNzeTHryFkGaK/GzUrBnnTaHFGb0ZqbiumJsDlsYpek5gc5BjAT07gttQG6vS2NgABRXNdJi08xLKfExDE3tuS9/QLgnOm0t2nJssma+iNBgOF9ON52mrk1Lrmw1WxmZ0m25+KAZkZ5AvFWbh6lvsVNS3dzrffQIk0nLiWV1ebi1NUZEUr/MuEzS47Q8XXannaN1R8MqT08JVEH8CdgsIg+JyIPAp8DzoRPLYCCz75R37IOE+kF9+H2ocznliUmblLaGsd5EkByobs+Un5+aj8UUaA7OwDGZhLG61Bv7whk050tssqYk3N+TmhNQvCGsIomIVxxKf/VmCrRg0CPA14BqoAb4mlLq0VAKZjAwabE5OFLRbl4aPzTE5qXyfd6V4fIXa26t/QSH0+EdHJd+9sFxnTF+SHuG/4Pl9dgdERQElp7rXbTp2MdQFd5IZv29OF5/nGZ7hIy6giDQEQRKqc+UUr90/X0eSqEMBi6HyhuwOzX7dnZyLFlJIUxr0XRGS/zmJqsAcmaHrr8QUNJQQotdM40lxiSGNPfP0NQ40hK0ye9Wm5OjlY0h66tHjJrv7VSwd51WRyJMpMamMihhEABO5eyXqTcCVhA9QUQuFpH9InJIRO7zs328iHwsIq0i8h2fbcUislNfSMgg+tGbLiaEcvTgsGvFf+yuqbX4tIisDNcdevNSQVpBSM1xIuI1itgbSWYmcLklL2t3S7a1uOYjHGETaVz6OM9n/b3qL4RMQYiIGfg1cAmaS+x1IuLrGnsG+B9gdSeHWayUKlRKzQqVnAaRQ12LjZJqLSJWBMYN6bJo4dlx6F0tjTdoQXATl2tVzPoRNoeNo7Xtk59jM3rfe8mX8TpvpuLKRprbwvfw9Ys10TUf4Xq01Z6EI+vDJs6Y9DEepV3WWEZNS03YZOkJoRxBzAEOKaWOKKXagBeBy/U7KKXKlVJbcCUBNBjYHCir9wTDjsxIICm29ydbASjbBaU6K2n+f0HK0ND0FUKO1B7xROmmx6WTGZcZ8j7TE60erzKHU3kVc4oY0kbA6EXtyyc2Q0V43t7jLfGMSh7lWe5vo4hQKojhgD6EsMS1LlAU8LaIbBOR23pVMoOIRG+yGB+q0UNjpXfxn0ETYPiM0PQVYvQPm7HpY0Pv7eVi/ND2e6P3OIsoRszV5pTc6NO29zHjMrzNTBETQxIAoVQQ/r6twVyZ+UqpGWgmqm+JyEK/nYjcJiJbRWRrRUVFT+Q0iADK61uorNdSOMSYhfxBISgMZG/T5h3cQUsJGTDukn437wBa3emSeq0+tiCMTQ+9ecnN2MFJmFzX7GRNM7VNEWgAENEKO8Wlasv2Vte97/u8SKNSRmE1a4WN6trqKGss63MZekooFUQJMEK3nAOUBtpYKVXq+l8OvIJmsvK337NKqVlKqVnZ2dlnIa5BONlTWuf5nJ+dRKyll5PjKQUH/62NIEArPKNP+NbP2H9mP8r1vpWTnOPJ/dMXJFgt5Ga1x4nsPhU+T6EuiYn3TrRYf1pLD97HWEwWr9Qb+6v397kMPSWUCmILUCAieSJiRSswtC6QhiKSKCLJ7s9otbB3hUxSg7DicCov76WJw0JgXjq1XZt7cDN2CSQN6v1++gCllNdDRm/C6Csm6sxMe0rrItdskjIM8s9vXz75GZze3edi6L2ZDtUc6jcZXkOmIJRSduAO4N/AXuBvSqndIrJSRFYCiMgQESkBvg18X0RKRCQFGAxsEpEdwGbgDaXUv/z3ZNDfOVrZ4PGGSY6zMCK9l6OY68vg4Lvty0OnwtBpvdtHH1LWWEZtq/bWbjVbyUvt+4SCo7OTvFJvnDgTwUFgw2d6F3va/1b7SLKPGJI4hBSrplTbHG1e3meRTIjcRDSUUm+ipQnXr3ta97kM/0n/6oD++ws2CIrdOvPSxKEpvZu51das2Z7db2xJ2VCwpPeOHwb2ndnn+TwmbQwxpt7L3BooZpMwfkgynx/X3DZ3l9YyMjNC05OIwLgvaW7NTWe0Oajdr8CMFWCx9pEIwriMcWwp06L2957ZG9Ko994ipIFyBgbd0dBq94rI7VXzklKw7w0tUytoD4NJV0AvpsLua2xOG4dr2yNyx2eErwyq/l4dKm/wJFiMSCyx2pyT2fVO3FgJB97q0yJD4zPGIy7fnZL6Ek+CxUjGUBAGYWXvqTrPbzQnPZ60hF58ozv+CVTqcvGPu1TzXOrHHKk5QptDi/5OjU1l8P9v78yD47rvw/75vj2wwGKxuA8CBAHeoMRTpCRaVuRIok3GthQ7jisnadIkMx5P42k6baZN6k47nbbTzCTTa+KJxznapvERx7FcST502JZlyTp4iCcoXiCJg7jva8/36x+/PcEFCBC72F3g95l53P29fb/dL7Hv7ff9vmdZQ95kqfd5qK/QTv6IrbhSaJnVCymvhx0fS44HO7VPYo3wuX20+JIGkw9GP1ji6MLAKAhD3lBKcakvGQHzwCZ/9t58/JbuDhen5Ui6HbpI6RztTDzfXb17zXIfFiP1O0s1FRYsTftg04Hk+MaPdLb1GrGnJllM4vLYZWxVQAUPM2AUhCFv9E8GGI/F0LudFtvry+8xY5kEpnQNnvjSxN+iq7QWOeOBcfpndVlyEcmreSnO7kYfjpjPaHAqwPB0AbQjvRfbj4GvUT+3o9ofEVqbwoNtFW14nDoTfTY8m8hlKVSMgjDkjfO9ydXDzgYfbmcWTsd485+QrumE25seC1/EpK4e2iva8bpykEy4QjwuR5piv9hXoDkRqThiOTDx2lvB6TVrMuSwHGkhr51jnUscnX+MgjDkhUA4yrWU3sZ7m7NkXrr+WtJkILFuYyVrl0SWKyJ2JC16KdVUkW9Sv7vO/ilCkcI2mwC6em/HM8ks+vHbcPP1NfnojpqOxPNbk7eYC8+tyefeD0ZBGPLCpTtTib4P9RUlNFRkIaO5/3y603HrR6Bqy2JHFxVdk10Eo9p843P72OzbfI8Za0dLVSlVsT4RoYjN1cECd1bHqdmme0jE6X5XN5DKMdWe6kRwga3sNMVfaBgFYVhzlFJc6E2WPd7fUrl6Z+v0AFx9OTmu3w2bM1ZnKUoujSSzfzuqO/LunE5FRNjbUpkYp5oOC562D0NNsgwGH7y0Jkl0D9Y+mHh+afRSwTqrjYIwrDm94/NpzumdDas0AYVm4eJ3kslw3lod0lpAP6KrYaFzOtVEUSg8sKkCZ4qzemAykGeJlokIdHxCm5xAJ9Fd/AfdbCiHbKvcRolDr5qnQ9N0T3Xn9PPuF6MgDGtO6h3mnqaK1Tmn7ah2MMZbS8aT4dYoQ3YtuDSaXD0UinN6IR6Xgx0piv58bxE1xnGVwoO/kkyimxuDyy/mNInOaTnTFP3F0cIsNWcUhGFNmQ1GuD40kxjvbVmlc/rGT7SDMU7HM+DNfeOctSIUDRWsc3oh+1K+y6uD04WdWb2Q8nrdcjbO6HW4+UZOP/KBmgcSmdU9Uz2J+lqFhFEQhjXlQt8kduzOrLmylNryVTinBy5A78nkuP3x9CYx64ArY1cSmdOVJZUF5ZxeSJPfQ51Pf5/hqCqOxLlU6jug9dHk+PbPc+q09pf42Vyhv0+FSvMzFQpGQRjWjEjUTjM97Nu8itXDZB9cSSnwW7czPSJlHaCU4sLIhcR4b+3egnJOL0RE0lYR53omsO0CLQO+GO1PQPXW5PiDF3UfiRzxYE3SWX157DJhu7CaLxkFYVgzrg7OMBtMlvXeUX+fzunApHYkpjqld39i3Til4/RM9zAR1ArV7XDnpe/DSuloqsDj0kmJk/NhukZm7jGjwLAs2PMMlFbpcTQCF78Nwdz8P1orWhPNnoLRINfGr91jxtpiFIRhTVBKcaY72RN4/+bKRImGFRGPMomXRnB5YO9nirYz3FKcHzmfeN5R3ZFoW1nIuBxW2iriTHcROavjuEph768mAx0CU3DpOzlpV2qJxd7avYnxueFzBdV8ySgIw5rQOz6fqNPjcsj9ZU4rpePU40t+sXTJhPjd3jpiIjCRCH0UJC1uvtDZ1+JP9qwen2dwqkhCXlPx1sCeX06uSif74OoPcxLZ1FHTkejpMR4Yp3u6cEJec6ogROS4iFwRkesi8ocZXt8tIm+LSFBE/mAlcw3Fxfs9yTvJVDPEirj5RrrTcMcxqGpbvXAFSOrqYUvFFvwlWax0m2N8Hhc7G5L1md4vxlUE6EzrrSlFHgcuQPfbWf+YEkdJWsjrueFzWf+M+yVnCkJEHMCXgRPAHuBzIrIwRm8M+GfAn97HXEORMDEXoms4acM92Hofd/z953VUSZzmh6D5UBakKzzmwnNcHr2cGO+t27vE0YXJoS3J7/jq4DQzweLowXwXmx/WJcLjdP0Uhi4vfvx9sq9uX1ozoZH5tW2Juhi5XEE8DFxXSnUppULAN4FnUw9QSg0ppU4CC13395xrKB7OdI8nVubttV6qvSu0pY/f1sv7ODXbYPvT2ROwwLgwcoGo0s782tJaWsozdeUtbBoqPDRXlgIQtRVni3UVIQI7j0Nla3Lf5Zdg6k5WP6bCXcHWymT01LmhwlhF5FJBNAM9KePe2L6szhWRz4vIKRE5NTw8fF+CGnLHbDDCpb5kPPyhla4eZke0g9COJV2V1+kKrdb6dJ+FoqG00NZD9YcKOrR1KQ5tSdZnOtc7UVyJc6lYDnjw08luhHYELvy9zrjOIgfqko2Mrk5cZSaU/wiwXF5lmc7q5Xp4lj1XKfVVpdRhpdThurq6ZQtnWBvOdI8nqrY2+j1sri5d/uTgNJz/VrIujtsbiy5ZfxFLcTpHO9NaiqbeVRYb2+rKqSnXq8VQxOZcT5GuIiAZ2RTvIRGa00oilL1S3Q3eBpq8TYCO+js7fDZr732/5FJB9AKpaZ8twHLXZauZaygQAuFoWt2lI23Vy78bjgT1BRivseRw6QvUUzzO2pUSsSNpDsoD9QewpHhXSiLC4S3JHuDv90wUR6+IxSirhgc/A1ZKzaaL39ah11niYP3BxPPO0c6894rI5dl3EtghIu0i4gaeA15Yg7mGAuH97uQPQm25m211yywyZ0fh0nfvDmetaMqRpIXBtfFrzIZ1fkeZsyyt81ixsqvRR0WpDuGcD0W5UAwd55aicjN0fDI9/DWL3ei2VGyhtrQW0DcM+V5F5ExBKKUiwBeBl4HLwLeUUpdE5Asi8gUAEWkUkV7gXwD/VkR6RaRisbm5ktWQfYKRKGdTTApH2pe5eojnOox1Jfft/Jh2TK9jonaU04OnE+P99ftxxu9UixiHJRxpS/qdztweJxIt4lUE6F4j255KjkeuZS1HQkQ43HA4Mb40con5yPyq3/d+yekZqJT6PvD9Bfu+kvJ8AG0+WtZcQ/FwoXcy4ZT0l7rYuZyyGkrplqGDKX162x6DTQcWn7NOuDJ+hamQduaXOEp4oOaBPEuUPfY0VfBu1xgzwQgzwQid/VPsS2kwVJRsPgLBKeh5T4/7z4G7THcxXCXt/naqPdWMBcYI22HOD5/nkaZHVv2+90PxGjgNBUswEuXkrWRZjSNt1VjLKatx+y3oPZUcbzoIbY/nQMLCImJHODWQ/H8fqD9QFGU1lovTYaVFNL13c4xwsa8iALY9CY0pOSq334aek4sfv0xEhIcaHkqML4xcIBDJTza6URCGrHP69nja6mHPpop7T+o9BTd/lhzX74YdH113Bfgy8cHYB8yEdUhjqbOUfbX77jGj+NjbXIm3RGfPTwcixdVQaDFEYNeJ9Jal11+DO6v3G2yr3EZliVaqoWgob74IoyAMWWUuFEkrrfDo1pp7F+W7cxauvZocV7XB7k+u21yHVMJ2OM33cLD+IC6HK48S5Qa30+Lh9mQjp5O3xglGijQvIhXLAQ/8MvhTLOVXfwiDq3OZWmJxuDHpizg/fD4RwLCWrP8r0LCmnLw1nha5tLvxHr6HgYvpWdL+5vT2j+ucSyOXEhe+1+Xlgdr143tYyN5mf1pE05nb62AVAckQbF+DHiuls62Hr67qbXdU7kiLaEo1Q64VRkEYssZUIMz5lMilo9tql/Y9DF2GD76XjP7wNcLez66rftJLEYgEODN0JjE+VH8oUdVzPeKwhKNbk6uIM93jzIfWwSoCdALdvud0bxIAZUPnd3WE030iIjzalOxw1znWyXhgfIkZ2ccoCEPWeOfGaFrW9JJ5D4Od0PmCvpBAX1j7/lEyU3UDcHrwdML56HP70ip6rld2N/rSsqvfuTmaZ4myiLsM9n8upSRHFC49vyolsdm3meZyXWVIKcW7A+9mQ9JlYxSEISsMTAbSehA/tq128byHwU64vEA57P+cvsA2CBOBibSaS0ebjq6LvId7YVnCh7bVJsbneyYZnQnmUaIsU1Kuz+V4j5JVKomFq4iuiS4GZgeyIemyMArCsGqUUrx+ZSgx3lrnpbVmkR/7gQsx5RAzK8WVQ0l55uPXKW/3v40dU5BN3ia2Va7vRMBUttV5aanSNblspfjp1eGC6qK2ajwVcODX0pXExe/cd5nwBm9D2vnxs76fJc6dXGMUhGHVXO6fpn9Sm0oclvDEzkWKJvae1s67Da4ceqd7uTl5MzH+0KYPFW3F1vtBRPjIrvpEBPPt0TluDK99hE5OWagklK1LcvTfXxnvR5sexSE6THh4bpjLY9nvSZEJoyAMqyIYifLW9WRzk0OtVVSWZXAy334brr2SHJfX6QtogymHqB3lrb63EuNdVbto8DbkUaL8UOcrSetd/cbV4eIvwbEQTwUc/PUUx7WCD76fngy6TPwl/rRCfu/ceWdNSnAYBWFYFe/dHEt0CysvcXKkfUG/B6Xg+o+g6/XkvopNcODXdfnuDcb5kfOMBrRj1mk581ZCoRA4urU20Xp2cj7M6dtrG6GzJpT49I2QL+Um4NqrujPdCs1qBxsO4nPrsPFgNMh7/e9lU9KMGAVhuG+GpgJpseyPba+lxJnSa9qOan9DT8qJXLUF9j+n6+tvMCaDk5wcSJZiONJ4hHL3xlpBpVLqdnB0WzLs9b2bY4zNhvIoUY5we2H/r+kcnzi3fw5XfrCiKrAuy8VjzY8lxp2jnQzODmZT0rswCsJwX0RtxSudg9ixu6DmqlI6mlKS4iJB3ewntfBe7Y5YnsP6bfizGEop3uh9g4itV1u1pbXsr9ufZ6nyz75mPw0VOrQ5Yite6xxcXw7rOPE8idSqxP3ndLfEyPKVYntFO60Vuv2pQvF6z+tE7dzlkhgFYbgvTt8eZ3hahyc6LeFYR0PS0To/AWf+BsZvJSc0H4IHPr1hMqQXcnX8Kj3TuouuIDzR8kRRNwPKFpYlPL2nHit27vRNzHOut8h7RiyG062rBKQW+Bu5Bmf/FgJTi89LQUR4vPnxREj0aGCUU4O5y7A2Z6hhxYzNhninK5ngdHRbDVXemGN6shfO/B/dSzpO+y/ownsboLZSJmbDs7zZ92ZivLdu74Z0TC9Gvc+T1jPiresjTM5nr0tbQWE5YPfHoTWZ28D0oL5mppeX3+Av8aflRpwZOsPQ3NASM+6fjXnFGu6bqK14+dIA0VjGdEOFh0OtsYt74AKc/UayT6/lgI5P6J4OGyiMMxWlFD/q/hHBqF5tlbvKeaRx4zqmF+Ph9uq0DOtXOwex7XVoagJ9LWz7Rdh1XHdLBAjOwPv/d9m5Entr96b1r/5x948T5stsklMFISLHReSKiFwXkT/M8LqIyP+MvX5eRA6lvHZLRC6IyFkRWfsqVYaMvH1jlIFYzoMlwrE9DVjYcPUVneMQP0ldpTrHIXU5vQE5N3yO3uleQJuWntry1Lqs1rpanA6LY3saEvcRPWNznFqPUU2pbDoI+1J8ctGIbrV748f3dF6LCE+2PpkwNY0FxnJiasqZghARB/Bl4ASwB/iciOxZcNgJYEds+zzw5wte/0Wl1AGl1GEMeef26Cwnb40lxo9tr6HOHYazX4e+ZMlqvLXw0G/p/r0bmJH5Ed7pfycxPthwMFFXx3A3Tf5SHm6rTozfvjFK30T+2m2uCdXtcOi3kvWbALrfhQvfgtDSyYP+Ej9Hm44mxpdGLiVWqtkilyuIh4HrSqkupVQI+Cbw7IJjngX+RmneASpFZH13pi9SZoMRfngxaSNtqy3jId8EnPor7XeIU7dLn/ClVRneZeMQioZ49fariZIIdWV1HGk4kmepCp9Ht9bQXJksw/GDC/2J5lPrFm+NvmZSGw+N3YRTfw0T3UtOfbD2QTaVb6LF18Jnd32WEkd2IwRzqSCagZ6UcW9s33KPUcArInJaRD6/2IeIyOdF5JSInBoeHs6C2IaFRG3FDy4OMBcrzex1C8fLriIXvpX0N4jofrwPfGrDlOtejLjfIV6a2Wk5OdZ6DIfluMdMg2UJx/c2JhLopgMRXr40sD5DX1NxeWDvZ7S/Lk5wRq/Ob/98UZOTiHCi/QSf3PrJRBJdNsmlgsjklVz4LS91zGNKqUNoM9TvicgvZPoQpdRXlVKHlVKH6+oWqQFkWBU/vTpEz5hWBJ7IJL8iP6X0zrvJTFC3V5fq3nJ0wzqjUzk1eCqt1tITLU9Q6alcYoYhlQqPi2N76hPjruFZ3rq+jsqCL4aIjvjb99lkIqlSOuv63Nd1+HgGShwlOavllUsF0QukGqFbgDvLPUYpFX8cAp5Hm6wMa8y5ngnO9UyCUtTPXObjwR9QE01ZqVW1weHf0bZUAzcnb6ZlS++v28+u6l15lKg42V7v46EtSTPlyVtjdN5ZXq5A0VOzTV9TqW1MJ3q0OXfgwopLdKyGXCqIk8AOEWkXETfwHPDCgmNeAH4zFs30KDCplOoXEa+I+ABExAt8FLiYQ1kNGegeneP1K8OURKbZPfwDjgTfZbM/lugmFmx9QpfN2GAF9xZjaG6I126/lhi3+Fo4uunoEjMMS/Hh7bW01ybrdb12eZA7691pHcdToeuVtT2WDIWNhHSk4MV/WHZi3WrJmYJQSkWALwIvA5eBbymlLonIF0TkC7HDvg90AdeBvwD+aWx/A/CmiJwD3gO+p5T6IYY1Y3AqwIvneqmfusC+/m/TogbYVl+OIMkopS0fMialGOOBcV7qeomwrRO8fG4fx7YcM9nSq8CyhBN7GxP5EVFb8cK5O4yspwZDS2FZ2uR08DfSgz5GrsHJv4A77+d8NSHryflz+PBhdeqUSZlYLSMzQb7/1hmaht6gPDSM22HxYHMFJS4nNB/WKwcTy59gOjTN89eeZyY8A2ib8Kd2fIpqT/U9ZhqWw+RcmG+c7E70r/aWOPjs4c2Zy8qvVyIh6PoJ9J1J3+9vhh0fS68Wu0JE5PRiqQTm9saQxvjEBKde+Ro7er9DeWgYpyXsbvJR4m/QdzI7njbKIYXZ8Cwv3ngxoRxclouPb/24UQ5ZxF/m4pcPNON26p+r2WCUb5/uZSqwTstxZMLphp0f0/0lUnMmJvvg9P/SJcTDgax/rFEQBk00wsSVN7j24p/iH+8EFA5L6GiuwrvzI/DQb6c7zQxMBid5/trzTAR1dIklFsfbjtPobcyzZOuPRr+HZw9swuXQJs3pQIRvn+plYm4dlgdfispW7cDeclSXsgFtZuo9pc1OkeJJlDMUA7YNAxeY+Omfce3nLxINxVqHCmzvOED5h78A7Y9v2CqsizEWGOO717/LVEg7C0WEY1uOsbliY2eP55KWqjI+sW8TDksricn5MH93soehqezfORc0DpfOOTr8uzqKME7NjqyX0jdX/UbFtmH4Mtx6i7GRfq4NThOvjRYuqaTt6DNUb9uXXxkLlL6ZPl6+9TKBSFyZOvhY28do87flV7ANQFutl0/u38T3zt8hHFXMhaL8/elePrlvE601ZfkWb23x1ugowuEPdDJde8ZUsVVhnNQbDTsKgxeh+13U3Ah94/P0js+jgIjlYbj2CB9+/GnqKzdeO9B7oZTi4shF3rzzZiKz12W5ONF+ghafMb+tJXcm5vnu2T6CYZ1hbInw4R21HGqtzFnSWEGj1H1HFC7lpDYriI1CeB76z0PvSQhOE7Ftrg/NMD4XJmq56PftY67uIM8e2bqxokOWSdgO82bvm1weS5ZjLnOWcaL9hOntkAc2VZby2cObef5MHzPBCLZSvHF1mMGpAE93NCQc2huGHClFoyDWOzNDutLq4EVdThhtu+0anmHWdjBQcZA7vr001Vbxmb1NeEvMKbGQeAJc3BkNuvjeibYTG7qndL6pLS/huYc389L5/kQJ+isD0wxPB/noAw00+Tde3/NsY34N1iORIAxegoHzMNWf2B21Fd1jc/TMCv2+wwyWdxC1SnhoSxUf3l6LZW3ApfkSRO0oZ4bOcGrwVFqxuJ1VO3li8xO4LBPum298Hhe/+lALP706zPlYq9Kx2RB/d7KHh7ZUcXRrDU7HBltNZBGjINYLdhRGb8DQJRi5nmzcg25uPjYb4oOZMm6VPMzwpu0ocVLisjje0cDOhuxXgSx2uqe6+Vnfz5hHLMqEAAANg0lEQVQMJvsjuywXjzU/Rkd1x8a0cxcoTofFUx0NNPo9vH5lmFDERik4dWucG0MzPL6zjq21XvOd3QdGQRQzkRCM34ThKzB6PWMM9HRIcSHYQKdsZ7q6MWGr3Frn5cnd9fg85i44lZH5Ed7rf49bU7fS9jd6G3mq9Sn8Jf78CGa4Jw9s8tNSVcarnYOJ6sPjc2FeOHuH1uoyHt9ZS73Pk2cpiwujIIoJpWB+XDcTGbsB47fTVgqJw1BMWFWci2zhgmohUpq8KDwuB7+ws5Y9TRXmjiqFkfkRTg2comuyK22/2+HmSOMR9tbuNXWVigB/qYtfOdTMhb5JfnZthFBERzl1j83xtXe62Vrn5eH2auOfWCZGQRQySkFgUndsm7itFUJgctHDoyUV9LjaeD+wiVuBWJhqLNnSYQn7N1fySHt1ohnLRidqR7k5eZOLoxe5M7OwEj3srt7No02PUubaYPH1RY6IsK+lku315bzTNcr53slETbuu4Vm6hmdpriplX4uf7XXlxkexBEZBFBKREMwMwPQATN3RiiE4veQUu6yGEU8bV+0mzk94CM6m57WIwM4GH0e31lDlNeGrSikG5wa5PnGd6+PXmYvM3XXMVv9WDjcepra0Ng8SGrJFmdvJk7sb2NdSyc9vjHJjaCbxWt/4PH3j83hcDnY3+djV4KPJ7zGr6gUYBZEPbBsCEzA3CrPDMDMIM8MwP3bv8r0OFyHfZgadm+iyG7gy6WB2Kt6zNznXYQkdTRUc3lK14RVDOBqmb6aP7ulubk/dZjp0t9IVEbb6t/JQw0NGMawzastLeGb/JkZmgpy6Nc6VgWns2HUWCEc52z3B2e4JykucbK8vp7WmjJaqUkqcZqVtMqlzRTQCwSmtCOYnYgphTPsQ5icy+g4yoRwu5kubGHHU0UcDt8MVDEyHF9UjlWUuHmz2s6epYkPmNCilmAnPMDw/zMDsAAOzAwzNDWGrzD19vS4ve2r20FHdYXIaNgjTgTCdd6a4eGeKqfnMFWEtERr9JTT5S2nye2jwe/CVONflCsNkUmcL24ZIQGclh+f0FpqNbTP6MTilzUKhu00XSxFVimBEMef0M+mqZcyqod+uoTdSTjCtFe3dJ7S3xMH2+nJ21PtoqSpdlyfxQiJ2hJnQDNOhacaD40wEJxgPjDMyP0IwunRFS7fDzTb/NrZXbae5vNk4nzcYPo+LR7bWcKStmp7xOa4NznB9eCbRbwLAVoo7EwHuTCQLAXpcDmrL3dSWl1BZ5qKyzE1lqQufx7lu/Rg5VRAichz4H2hX6V8qpf54wesSe/2XgDngnyilzixnbtYYvaHv6KOh5BYJQjSofQKRQMoWvKcJSClFVClsBbatiNiKqK2I2jaR2DgSVcyLhxlHJZOWnzHxM0YV865q7KgTokt+BCJQ7/OwpaaMLTVlbPKXFm2Sm1KKiB0hbIcJ22FC0RAhO0QwGiQUDRGIBAhEA8yH55mLzDEbnmUuPJfRd7AU1Z5qWn2ttFa00uRtwmEZ88FGx7KELTVettR4edKup29inlujs3SPzTE0dfdNRiAcpTdWu2whZW4H5R4nXreTMreDMreTUrdFidOBx+WgxGlR4rRwOy1cDgunQ3BZVsFftzlTECLiAL4MHAN6gZMi8oJSqjPlsBPAjtj2CPDnwCPLnJsVXnn3G4xO3NQDlbTiqxR7vlIkMmlVfIv9Y6Nf08dAqiEjbHkIO0oJW6WEHWWEYlvY6cVOzcJVM8AMhHoyyuhyCP4yNzVlLqq8bqq9btxOB1Gga1ZvyzUVpv6/0vanzI8fkzhWJfcrVOJYGzv2N9P7bGz9t0BhK/3cVjZRFU08Rm39PKIiROzIoqaf1eB2uKnx1NBQ1kBjeSONZY0mEsmwJJYlbK4uY3O1Pk/mQhHuTAQYmAzQPznP0HQwETKbiblQlLlQFFhZPwaHJTgdgtMSHJaFQ/Q+yxIsERwiiGiTl2WBoMcigkDiNUH3zDjYWnWvj1wRuVxBPAxcV0p1AYjIN4FngdQf+WeBv1H6F+cdEakUkSagbRlzs8IHc+P0B0eXfbwSB7Y4iOLEtlxExYktLqKWi2ji0U3UcgEL7w7mQM1BZOSu9xUBt9PC43TgcVl4XA5K3Q7K3A7cDgsRIQqMhGFk4q7pGw4Rwev04nP78Jf4qfJU4S/xU+OpocJtcjwMq6PMrR3W2+u1X0opxVQgwshMkPHZEONzYSbmQkzOh5kNRhNO75USjVkYtFq5h+ngHijg4Kre4W5yqSCagdRb4l70KuFexzQvcy4AIvJ54PMAra2tKxYy5KlhNjCFLRZKLBRaAdjiSCgDvTmxxcHdP/oL5dEa3RW7E0jdnJaVuFtwOfRS0+2ILztlw/6oOS0nLsuV2NwON26HmxJHCR6nB4/DQ6mzlFJXKWXOMrwuL16X1/gODGuGiOAvdeEvdUFd+mu2rZgNRZgJRpgNRpkPRZkLRZgPRwmEbYKRKMGITSi2haPa3LzUiuS+ZMzqu2lyqSAyybtQzS52zHLm6p1KfRX4KugoppUICPDYvk+za24iuXRLkcqK/djr57GlniVYsWMdsWWgZQnOlKWhrPCrWgvFsKhMsvhx8ed6OZuyP/43wUq+Jslx4lGstM0hDiyxcFpOHOLAYTlwyvqMDDFsHCxL8HlcKy5bo5RePcR9k9Go9l9GbYWt9Ba1tfnajvk1VewRYmZtYiZwVE7K5uRSQfQCqf0XW4CF6aqLHeNextyscLR1dy7e1mAwGJZEJOZ/KOB4iVyu0U8CO0SkXUTcwHPACwuOeQH4TdE8CkwqpfqXOddgMBgMOSRnKwilVEREvgi8jA5V/Wul1CUR+ULs9a8A30eHuF5Hh7n+9lJzcyWrwWAwGO7GZFIbDAbDBmapTGoTBmIwGAyGjBgFYTAYDIaMGAVhMBgMhowYBWEwGAyGjKwrJ7WIDAO373N6LXB3DYz8Y+RaGUaulWHkWhnrUa4tSqm6TC+sKwWxGkTk1GKe/Hxi5FoZRq6VYeRaGRtNLmNiMhgMBkNGjIIwGAwGQ0aMgkjy1XwLsAhGrpVh5FoZRq6VsaHkMj4Ig8FgMGTErCAMBoPBkBGjIAwGg8GQEaMgMiAifyAiSkRq8y0LgIj8RxE5LyJnReQVEdmUb5kARORPROSDmGzPi0hlvmUCEJFfFZFLImKLSF5DEkXkuIhcEZHrIvKH+ZQlFRH5axEZEpGL+ZYlFRHZLCI/EZHLse/w9/MtE4CIeETkPRE5F5PrP+Rbpjgi4hCR90XkpWy/t1EQCxCRzcAxoDvfsqTwJ0qpfUqpA8BLwL/Lt0AxXgUeVErtA64Cf5RneeJcBD4NvJFPIUTEAXwZOAHsAT4nInvyKVMK/xs4nm8hMhAB/qVSqgN4FPi9AvmbBYEnlVL7gQPA8VgPm0Lg94HLuXhjoyDu5r8B/4pFWpzmA6XUVMrQS4HIppR6RSkViQ3fQXf+yztKqctKqSv5lgN4GLiulOpSSoWAbwLP5lkmAJRSbwBj+ZZjIUqpfqXUmdjzafQPX3N+pQKlmYkNXbEt79ehiLQAHwf+MhfvbxRECiLyDNCnlDqXb1kWIiL/WUR6gF+ncFYQqfwO8IN8C1FgNAM9KeNeCuDHrlgQkTbgIPBufiXRxEw5Z4Eh4FWlVCHI9d/RN7R2Lt48lz2pCxIReQ1ozPDSl4B/A3x0bSXSLCWXUur/KaW+BHxJRP4I+CLw7wtBrtgxX0KbBr62FjItV64CQDLsy/tdZzEgIuXAPwD/fMEKOm8opaLAgZiv7XkReVAplTcfjoh8AhhSSp0WkY/k4jM2nIJQSj2dab+I7AXagXMiAtpcckZEHlZKDeRLrgx8Hfgea6Qg7iWXiPwW8AngKbWGSTUr+Hvlk15gc8q4BbiTJ1mKBhFxoZXD15RS38m3PAtRSk2IyOtoH04+nfyPAc+IyC8BHqBCRP5WKfUb2foAY2KKoZS6oJSqV0q1KaXa0Bf3obVQDvdCRHakDJ8BPsiXLKmIyHHgXwPPKKXm8i1PAXIS2CEi7SLiBp4DXsizTAWN6LuzvwIuK6X+a77liSMidfEoPREpBZ4mz9ehUuqPlFItsd+r54AfZ1M5gFEQxcIfi8hFETmPNoEVROgf8GeAD3g1FoL7lXwLBCAinxKRXuAo8D0ReTkfcsQc+F8EXkY7W7+llLqUD1kWIiLfAN4GdolIr4j8br5livEY8I+BJ2Pn1NnYHXK+aQJ+ErsGT6J9EFkPKy00TKkNg8FgMGTErCAMBoPBkBGjIAwGg8GQEaMgDAaDwZARoyAMBoPBkBGjIAwGg8GQEaMgDAaDwZARoyAMBoPBkBGjIAyGHCEiR2K9Mjwi4o31EXgw33IZDMvFJMoZDDlERP4Tuk5OKdCrlPoveRbJYFg2RkEYDDkkVoPpJBAAPhSrCGowFAXGxGQw5JZqoBxds8qTZ1kMhhVhVhAGQw4RkRfQneTagSal1BfzLJLBsGw2XD8Ig2GtEJHfBCJKqa/H+lP/XESeVEr9ON+yGQzLwawgDAaDwZAR44MwGAwGQ0aMgjAYDAZDRoyCMBgMBkNGjIIwGAwGQ0aMgjAYDAZDRoyCMBgMBkNGjIIwGAwGQ0b+P+nNAS3LrWZNAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(support, init_untwist_dist, lw=3, alpha=.5, label = f\"t=0 untwisted distribution\")\n", + "plt.plot(support, init_twist_dist, lw=3, alpha=.5, label=f\"t=0 twisted distribution\")\n", + "plt.plot(support, final_dist, lw=3, alpha=.5, label=f\"t=9 (last) final target distribution\")\n", + "plt.xlabel(f\"x\")\n", + "plt.ylabel(f\"density\")\n", + "plt.title(f\"initial and final target distributions\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "cool" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "NOTES:\n", + "- in both of these examples, the twisting functions that control the variance and mean of the propagations are _independent_ of the position of the latent variable, which is perhaps an unnecessary constraint\n", + "- to fix this ^, I will write more expressive functions (while maintaining the constraint that the covariance must be positive-definite)\n", + "- I need to fix some of the loggings since they are perhaps a little overly verbose\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/melange/smc.py b/melange/smc.py index 43218df..a8d846e 100644 --- a/melange/smc.py +++ b/melange/smc.py @@ -7,6 +7,7 @@ from jax.scipy.special import logsumexp import numpy as np from melange.reporters import BaseSMCReporter +import tqdm import jax from jax.config import config; config.update("jax_enable_x64", True) @@ -180,3 +181,26 @@ def scanner(none, t): out_works = jnp.vstack((init_logWs, logW_matrix)) return cond(aggregate_works, lambda x: jnp.cumsum(x, axis=0), lambda x: x, out_works) + +def bootstrap_logZ(incremental_work_array, num_bootstraps=1000): + """ + bootstrap resample the logZ estimate + + arguments + incremental_work_array : np.array(T,N) + work array (unaggregated) + num_bootstraps : int, default 1000 + number of bootstraps to conduct + + returns + bootstraps : np.array(num_bootstraps) + logZ estimates + """ + T, N = incremental_work_array.shape + final_works = np.cumsum(incremental_work_array, axis=0)[-1] + bootstraps = np.zeros(num_bootstraps) + for i in tqdm.trange(num_bootstraps): + _works = np.random.choice(final_works, size=N, replace=True) + logZ = logsumexp(_works) - np.log(N) + bootstraps[i] = logZ + return bootstraps diff --git a/melange/tests/test_csmc.py b/melange/tests/test_csmc.py index 9195243..b0d8b82 100644 --- a/melange/tests/test_csmc.py +++ b/melange/tests/test_csmc.py @@ -124,6 +124,37 @@ def test_null_do_param_twist(dim=3, num_twists=5, scale=1e3): assert np.isclose(log_normalizer, 0.), f"{log_normalizer}" +def test_log_psi_twist(dim=3, num_twists=5): + """ + test the log_psi_twist function + """ + x = np.random.randn(dim) + As = np.random.randn(num_twists, dim) + bs = np.random.randn(num_twists, dim) + + out = log_psi_twist(x, As, bs) + assert type(out) == np.float64 + + #manual + outs_manual = np.array([np.dot(x*A, x) + np.dot(x,b) for A,b in zip(As, bs)]) + assert np.isclose(-outs_manual.sum(), out) + +def test_vlog_psi_twist_single(num_particles = 10, dim=3): + """ + test the vlog_psi_twist_single function + """ + xs = jnp.array(np.random.randn(num_particles, dim)) + A = jnp.array(np.random.randn(dim)) + b = jnp.array(np.random.randn(dim)) + + out_twists = vlog_psi_twist_single(xs, A, b) + + #manual + manual_outs = np.array([log_psi_twist(x, A[jnp.newaxis, ...], b[jnp.newaxis]) for x in xs]) + + #make the assertions of isclose + assert np.allclose(out_twists, manual_outs) + def test_uncontrolled_csmc(dim=1): """ test a vanilla implementation of uncontrolled sequential monte carlo of a static model with an Euler Maruyama forward/backward kernel.