diff --git a/README.md b/README.md index ca3f7bc..23db213 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,80 @@ Copyright (c) 2023, University of Liverpool. - pandas - scipy - scikit-learn (for examples) + - [bridgestan](https://github.com/roualdes/bridgestan) (for models with continuous parameters specified with Stan) + - cmdstanpy (for models with continuous parameters specified with Stan) + +## Continuous Parameters and Stan +For models with continuous parameters, the continuous model conditional on discrete +parameters can be specified using a model written in Stan. +Log density and gradient information are extracted from Stan models using +[bridgestan](https://github.com/roualdes/bridgestan). + + + +## Installation + +Install with pip: +```bash +pip install pip install git+https://github.com/DiscreteVariablesTaskForce/DiscreteSamplingFramework.git +``` + +or + +``` +git clone --recurse-submodules https://github.com/DiscreteVariablesTaskForce/DiscreteSamplingFramework.git +pip install -e DiscreteSamplingFramework +``` + +### bridgestan + +[bridgestan](https://github.com/roualdes/bridgestan) must be "installed". +It is set as a submodule of this repo and so, if you installed by cloning the repo, you can +set an environment variable to point to the bridgestan repo: +```bash +export BRIDGESTAN=/bridgestan +``` + + +Alternatively, you can clone bridgestan elsewhere and set an environment variable: +```bash +git clone https://github.com/roualdes/bridgestan.git +export BRIDGESTAN= +``` + +If this environment variable is not set, discretesampling will attempt to set it to +the bridgestan submodule. + + +### Cmdstan +[cmdstan](https://github.com/stan-dev/cmdstan) is also required to query stan models with bridgestan. + +#### Installation with cmdstanpy +The simplest way to install cmdstan is using cmdstanpy: +```python +import cmdstanpy +cmdstanpy.install_cmdstan() +``` +which will download and build the latest version of cmdstan to +your home directory in `~/.cmdstan` and +sets an environment variable `CMDSTAN` to that location. + +#### Manual installation +Alternatively, cmdstan can be installed manually following instructions for [installing with conda](https://mc-stan.org/docs/cmdstan-guide/cmdstan-installation.html#conda-install) or [building from source](https://mc-stan.org/docs/cmdstan-guide/cmdstan-installation.html#installation-from-github). + +The `CMDSTAN` environment variable should be then set: +```bash +export CMDSTAN=path/to/cmdstan +``` +or with cmdstanpy: +```python +import cmdstanpy +cmdstanpy.set_cmdstan_path('path/to/cmdstan') +``` + +#### Windows +On Windows there's an additional step in the cmdstan installation to link the TBB libraries. However, this might not work correctly. In this case you may run into an error when using bridgestan, where it will say that it is unable to find the model .so file or one of its dependencies. The way to fix this is to copy the tbb.dll file from `cmdstan\stan\lib\stan_math\lib\tbb` to the folder containing the stan model. + ### Cloning and installing from github diff --git a/discretesampling/base/algorithms/continuous/NUTS.py b/discretesampling/base/algorithms/continuous/NUTS.py new file mode 100644 index 0000000..0b4272b --- /dev/null +++ b/discretesampling/base/algorithms/continuous/NUTS.py @@ -0,0 +1,306 @@ +import copy +import numpy as np +from scipy.stats import multivariate_normal +from discretesampling.base.algorithms.continuous.base import ContinuousSampler + +MAX_TREEDEPTH = 10 + + +class NUTS(ContinuousSampler): + # We need to adapt the NUTS parameters before we sample from a specific model as defined by the discrete parameters. + # Each time we jump to a new model, initialise NUTS to adapt the mass matrix and step size parameters and then store these + # for later use. + def __init__(self, do_warmup, model, data_function, rng, adapt_delta=0.9, warmup_iters=100): + self.NUTS_params = {} + self.current_stepsize = None + self.do_warmup = do_warmup + self.stan_model = model + self.data_function = data_function + self.rng = rng + self.delta = adapt_delta + self.warmup_iters = warmup_iters + + # Perform a single leapfrog step + def Leapfrog(self, current_data, theta, r, epsilon): + [L, dtheta] = self.stan_model.eval(current_data, theta) + r1 = r + epsilon*np.array(dtheta)/2.0 + theta1 = theta + epsilon*r + [L1, dtheta_list] = self.stan_model.eval(current_data, theta1) + r1 = r1 + epsilon*np.array(dtheta_list)/2.0 + + return theta1, r1, L1, L + + # Recursive part of NUTS algorithm which decides direction of tree growth + def BuildTree(self, current_data, theta, r, u, v, j, epsilon, treedepth, init_energy, init_potential, M): + treedepth += 1 + if treedepth > MAX_TREEDEPTH: + print("max tree depth exceeded") + return 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 + Delta_max = 1000 + n1 = 0 + s1 = 0 + if j == 0: + # take one leapfrog step + [theta1, r1, L1] = self.Leapfrog(current_data, theta, r, v*epsilon)[0:3] + energy1 = L1 - 0.5 * np.dot(np.linalg.solve(M, r1), r1) - init_potential + if u <= np.exp(energy1): + n1 = 1 + if u == 0: + # stop log(0) error + s1 = 1 + else: + if energy1 > (np.log(u) - Delta_max): + s1 = 1 + theta_n = theta1 + r_n = r1 + theta_p = theta1 + r_p = r1 + alpha1 = np.exp(energy1 - init_energy) + if alpha1 > 1: + alpha1 = 1 + n_alpha1 = 1 + else: + # build the left and right subtrees using recursion + [theta_n, r_n, theta_p, r_p, theta1, n1, s1, alpha1, n_alpha1, r1, depth_exceeded] = \ + self.BuildTree(current_data, theta, r, u, v, j-1, epsilon, + treedepth, init_energy, init_potential, M) + if depth_exceeded == 1: + return 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 + if s1 == 1: + if v < 0: + [theta_n, r_n, _, _, theta2, n2, s2, alpha2, n_alpha2, r2, depth_exceeded] = \ + self.BuildTree(current_data, theta_n, r_n, u, v, j-1, epsilon, + treedepth, init_energy, init_potential, M) + else: + [_, _, theta_p, r_p, theta2, n2, s2, alpha2, n_alpha2, r2, depth_exceeded] = \ + self.BuildTree(current_data, theta_p, r_p, u, v, j-1, epsilon, + treedepth, init_energy, init_potential, M) + if depth_exceeded == 1: + return 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 + if n1 != 0 and n2 != 0: + u1 = self.rng.uniform(0, 1) + if u1 < n2/(n1+n2): + theta1 = theta2 + r1 = r2 + alpha1 = alpha1 + alpha2 + n_alpha1 = n_alpha1 + n_alpha2 + arg = theta_p - theta_n + I1 = 0 + if np.dot(arg, r_n) >= 0: + I1 = 1 + I2 = 0 + if np.dot(arg, r_p) >= 0: + I2 = 1 + s1 = s2*I1*I2 + n1 = n1 + n2 + + return theta_n, r_n, theta_p, r_p, theta1, n1, s1, alpha1, n_alpha1, r1, 0 + + def NUTS(self, current_continuous, current_discrete, M, epsilon): + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + + # calculate likelihood for starting state + L = self.stan_model.eval(current_data, current_continuous[0:param_length])[0] + + # randomly sample momenta (mass matrix not currently implemented) + r0 = self.rng.randomMvNormal(np.zeros([param_length]), np.identity(param_length)) + + # calculate energy + init_energy = -0.5 * np.dot(np.linalg.solve(M, r0), r0) + # more numerically-stable if energy is normalised relative to initial + # logp value (which then needs to be passed into BuildTree()) + # note there should also be a "- np.log(np.linalg.det(M))" term here but + + # it's constant and so doesn't affect the algorithm + u = self.rng.uniform(0, np.exp(init_energy)) + theta_n = copy.deepcopy(current_continuous[0:param_length]) + theta_p = copy.deepcopy(current_continuous[0:param_length]) + p_params = copy.deepcopy(current_continuous) + r_n = r0 + r_p = r0 + j = 0 + s = 1 + n = 1 + treedepth = 0 + + # start building tree + while s == 1: + v_j = self.rng.uniform(-1, 1) + if v_j < 0: + [theta_n, r_n, _, _, theta1, n1, s1, alpha, n_alpha, r1, depth_exceeded] = \ + self.BuildTree(current_data, theta_n, r_n, u, v_j, j, epsilon, treedepth, init_energy, L, M) + else: + [_, _, theta_p, r_p, theta1, n1, s1, alpha, n_alpha, r1, depth_exceeded] = \ + self.BuildTree(current_data, theta_p, r_p, u, v_j, j, epsilon, treedepth, init_energy, L, M) + if depth_exceeded == 1: + # max tree depth exceeded, restart from beginning + r0 = self.rng.randomMvNormal(np.zeros([param_length]), np.identity(param_length)) + init_energy = -0.5 * np.dot(np.linalg.solve(M, r0), r0) + u = self.rng.uniform(0, np.exp(init_energy)) + theta_n = copy.deepcopy(current_continuous[0:param_length]) + theta_p = copy.deepcopy(current_continuous[0:param_length]) + p_params = copy.deepcopy(current_continuous) + r_n = r0 + r_p = r0 + j = 0 + s = 1 + n = 1 + treedepth = 0 + else: + if s1 == 1: + u1 = self.rng.uniform(0, 1) + if n1/n > u1: + p_params[0:param_length] = theta1 + n += n1 + u1 = 0 + if np.dot((theta_p - theta_n), r_n) >= 0: + u1 = 1 + u2 = 0 + if np.dot((theta_p - theta_n), r_p) >= 0: + u2 = 1 + s = s1*u1*u2 + j += 1 + + return p_params, r0, r1, alpha, n_alpha + + def FindReasonableEpsilon(self, current_continuous, current_discrete): + epsilon = 1 + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + + # randomly sample momenta + r = self.rng.randomMvNormal(np.zeros([param_length]), np.identity(param_length)) + + [_, _, L1, L] = self.Leapfrog(current_data, current_continuous[0:param_length], r, epsilon) + + # work out whether we need to increase or decrease step size + u = 0 + if (L1 - L) > np.log(0.5): + u = 1 + a = 2*u - 1 + + # iteratively adjust step size until the acceptance probability is approximately 0.5 + while a*(L1 - L) > -a*np.log(2): + epsilon = (2**a)*epsilon + theta1 = self.Leapfrog(current_data, current_continuous[0:param_length], r, epsilon)[0] + L1 = self.stan_model.eval(current_data, theta1)[0] + + return epsilon + + def adapt(self, current_continuous, current_discrete, log_epsilon, log_epsilon_bar, M, H_bar, n): + mu = np.log(10) + log_epsilon + gamma = 0.05 + t0 = 10 + kappa = 0.75 + [proposed_continuous, r0, r1, alpha, n_alpha] = self.NUTS(current_continuous, current_discrete, M, + np.exp(log_epsilon)) + H_bar = (1 - 1 / (n + t0)) * H_bar + (self.delta - alpha / n_alpha) / (n + t0) + log_epsilon = mu - np.sqrt(n) * H_bar / gamma + n_scaled = n ** (-kappa) + log_epsilon_bar = n_scaled * log_epsilon + (1 - n_scaled) * log_epsilon_bar + + return proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar + + def warmup(self, current_continuous, current_discrete, log_epsilon, M): + # Step 2: Adapt step size using dual averaging + H_bar = 0 + log_epsilon_bar = 0 + proposed_continuous = copy.deepcopy(current_continuous) + for n in range(1, self.warmup_iters): + [proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar] = self.adapt(proposed_continuous, + current_discrete, log_epsilon, + log_epsilon_bar, M, H_bar, n) + + return proposed_continuous, r0, r1, log_epsilon_bar + + def continual_adaptation(self, current_continuous, current_discrete, log_epsilon, log_epsilon_bar, M, H_bar, n): + [proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar] = self.adapt(current_continuous, current_discrete, + log_epsilon, log_epsilon_bar, M, + H_bar, n) + + return proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar + + def init_NUTS(self, current_continuous, current_discrete): + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + M = np.identity(param_length) # currently no mass matrix adaptation (assume identity matrix) + + # Step 1: Get initial estimate for optimal step size + if self.current_stepsize is None: + log_epsilon = np.log(self.FindReasonableEpsilon(current_continuous, current_discrete)) + else: + log_epsilon = np.log(self.current_stepsize) + + if self.do_warmup: + # adapt parameters in one go so no need to save the dual averaging variables + [proposed_continuous, r0, r1, log_epsilon_bar] = self.warmup(current_continuous, current_discrete, log_epsilon, M) + if hasattr(current_discrete.value, "__iter__"): + self.NUTS_params[','.join(current_discrete.value)] = (M, np.exp(log_epsilon_bar)) + else: + self.NUTS_params[str(current_discrete.value)] = (M, np.exp(log_epsilon_bar)) + else: + # do one iteration of adaptation + H_bar = 0 + log_epsilon_bar = 0 + n = 1 + [proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar] = \ + self.continual_adaptation(current_continuous, current_discrete, log_epsilon, log_epsilon_bar, M, H_bar, n) + if hasattr(current_discrete.value, "__iter__"): + self.NUTS_params[','.join(current_discrete.value)] = (M, np.exp(log_epsilon_bar), log_epsilon, H_bar, n + 1) + else: + self.NUTS_params[str(current_discrete.value)] = (M, np.exp(log_epsilon_bar), log_epsilon, H_bar, n + 1) + + # save step size for initialising other parameter spaces + self.current_stepsize = np.exp(log_epsilon_bar) + + return proposed_continuous, r0, r1 + + def sample(self, current_continuous, current_discrete): + if str(current_discrete.value) in self.NUTS_params.keys(): + if self.do_warmup: + # run for a set amount of warmup iterations in order to optimise NUTS + if hasattr(current_discrete.value, "__iter__"): + [M, epsilon] = self.NUTS_params[','.join(current_discrete.value)] + else: + [M, epsilon] = self.NUTS_params[str(current_discrete.value)] + [proposed_continuous, r0, r1] = self.NUTS(current_continuous, current_discrete, M, epsilon)[0:3] + else: + # continually optimise NUTS (e.g. for SMC) + if hasattr(current_discrete.value, "__iter__"): + [M, epsilon, log_epsilon, H_bar, n] = self.NUTS_params[','.join(current_discrete.value)] + [proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar] = \ + self.continual_adaptation(current_continuous, current_discrete, + log_epsilon, np.log(epsilon), M, H_bar, n) + self.NUTS_params[','.join(current_discrete.value)] = (M, np.exp(log_epsilon_bar), + log_epsilon, H_bar, n + 1) + else: + [M, epsilon, log_epsilon, H_bar, n] = self.NUTS_params[str(current_discrete.value)] + [proposed_continuous, r0, r1, log_epsilon, log_epsilon_bar, H_bar] = \ + self.continual_adaptation(current_continuous, current_discrete, + log_epsilon, np.log(epsilon), M, H_bar, n) + self.NUTS_params[str(current_discrete.value)] = (M, np.exp(log_epsilon_bar), + log_epsilon, H_bar, n + 1) + else: + [proposed_continuous, r0, r1] = self.init_NUTS(current_continuous, current_discrete) + + return proposed_continuous, r0, r1 + + def eval(self, current_continuous, current_discrete, proposed_continuous, current_r, proposed_r): + # pull mass matrix from the NUTS parameters (needed to calculate kinetic energy) + if hasattr(current_discrete.value, "__iter__"): + M = self.NUTS_params[','.join(current_discrete.value)][0] + else: + M = self.NUTS_params[str(current_discrete.value)][0] + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + init_L = self.stan_model.eval(current_data, current_continuous[0:param_length])[0] + L = self.stan_model.eval(current_data, proposed_continuous[0:param_length])[0] + + log_acceptance_ratio = L - 0.5 * np.dot(np.linalg.solve(M, proposed_r), proposed_r) \ + - init_L + 0.5 * np.dot(np.linalg.solve(M, current_r), current_r) + + if log_acceptance_ratio > 0: + log_acceptance_ratio = 0 + + return log_acceptance_ratio diff --git a/discretesampling/base/algorithms/continuous/RandomWalk.py b/discretesampling/base/algorithms/continuous/RandomWalk.py new file mode 100644 index 0000000..335f194 --- /dev/null +++ b/discretesampling/base/algorithms/continuous/RandomWalk.py @@ -0,0 +1,40 @@ +import copy +import numpy as np +from scipy.stats import multivariate_normal +from discretesampling.base.algorithms.continuous.base import ContinuousSampler + + +class RandomWalk(ContinuousSampler): + def __init__(self, model, data_function, rng): + self.stan_model = model + self.data_function = data_function + self.rng = rng + + def sample(self, current_continuous, current_discrete): + proposed_continuous = copy.deepcopy(current_continuous) + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + mu = [0 for i in range(param_length)] + sigma = np.identity(param_length) * 1 + proposed_continuous[0:param_length] = current_continuous[0:param_length] + self.rng.randomMvNormal(mu, sigma) + + return proposed_continuous, 0, 0 # zeros added for consistency with NUTS method + + def eval(self, current_continuous, current_discrete, proposed_continuous, r0, r1): # r0 and r1 added for consistency + current_data = self.data_function(current_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + sigma = np.identity(param_length) * 1 + + current_target = self.stan_model.eval(current_data, current_continuous[0:param_length])[0] + proposed_target = self.stan_model.eval(current_data, proposed_continuous[0:param_length])[0] + forward_logprob = multivariate_normal.logpdf(proposed_continuous[0:param_length], + mean=current_continuous[0:param_length], cov=sigma) + reverse_logprob = multivariate_normal.logpdf(current_continuous[0:param_length], + mean=proposed_continuous[0:param_length], cov=sigma) + # Discrete part of target p(discrete_variables) cancels + log_acceptance_ratio = (proposed_target - current_target + + reverse_logprob - forward_logprob) + if log_acceptance_ratio > 0: + log_acceptance_ratio = 0 + + return log_acceptance_ratio diff --git a/discretesampling/base/algorithms/continuous/base.py b/discretesampling/base/algorithms/continuous/base.py new file mode 100644 index 0000000..97f406d --- /dev/null +++ b/discretesampling/base/algorithms/continuous/base.py @@ -0,0 +1,60 @@ +from abc import ABC, abstractmethod +from typing import Callable +import numpy as np +from discretesampling.base.types import DiscreteVariable +from discretesampling.base.random import RNG +from discretesampling.base.stan_model import StanModel + + +class ContinuousSampler(ABC): + """Abstract base class for continuous samplers + + Parameters + ---------- + model : StanModel + Stan model desribing continuous posterior + data_function : Callable + Function which generates Stan model data from a :class:`DiscreteVariable` + rng : RNG, optional + RNG for random number generation, by default RNG() + """ + + def __init__(self, model: StanModel, data_function: Callable, rng: RNG = RNG(), **kwargs): + """Constructor method + """ + self.model = model + self.data_function = data_function + self.rng = rng + + @abstractmethod + def sample( + self, + current_continuous: np.ndarray, + current_discrete: DiscreteVariable, + **kwargs + ) -> np.ndarray: + """Sample vector of parameters in continuous space + + Parameters + ---------- + current_continuous : np.ndarray + Current array of parameters in continuous space + current_discrete : DiscreteVariable + Current value in discrete space + + Returns + ------- + np.ndarray + Sampled parameters in continuous space + """ + pass + + @abstractmethod + def eval( + self, + current_continuous: np.ndarray, + current_discrete: DiscreteVariable, + proposed_continuous: np.ndarray, + **kwargs + ) -> np.ndarray: + pass diff --git a/discretesampling/base/algorithms/continuous_proposals.py b/discretesampling/base/algorithms/continuous_proposals.py new file mode 100644 index 0000000..6861716 --- /dev/null +++ b/discretesampling/base/algorithms/continuous_proposals.py @@ -0,0 +1,281 @@ +import copy +import numpy as np +from scipy.special import logsumexp +from ...base.random import RNG +from ...base.util import u2c, c2u_array + + +def sample_offsets(grid_size, min_vals, max_vals, rng=RNG()): + # input: + # grid_size - number of points across which to sample each new parameter + # min_vals - list of minimum values of each constrained parameter in order that they appear in + # current_continuous + # max_vals - list of maximum values of each constrained parameter in order that they appear in + # current_continuous + + num_params = len(min_vals) + offsets = [] + for n in range(0, num_params): + step = (max_vals[n] - min_vals[n]) / (grid_size + 1) + offsets.append(step * rng.random()) # random offset so that any point in the parameter space can be sampled + + return offsets + + +def grid_search_logprobs(data_function, stan_model, grid_size, min_vals, max_vals, offsets, inds, current_continuous, + proposed_discrete): + # input: + # data_function - function returning data passed to Stan model for given discrete parameters + # stan_model - initialised stan_model object + # grid_size - number of points across which to sample each new parameter + # min_vals - list of minimum values of each constrained parameter in order that they appear in + # current_continuous + # max_vals - list of maximum values of each constrained parameter in order that they appear in + # current_continuous + # offsets - randomly sampled offsets used to define the grids + # inds - indices of parameters to sample (starting with 0) + # current_continuous - current unconstrained parameters (including those with unchanging dimensions) + # proposed_discrete - proposed discrete parameter (assumed to be 1 parameter controlling birth / death moves) + + dim = proposed_discrete.value + + proposed_data = data_function(proposed_discrete) + + param_length = stan_model.num_unconstrained_parameters(proposed_data) + + num_params = len(min_vals) + + total_params = len(current_continuous) + num_static_params = total_params - (dim - 1) * num_params + if num_static_params < 0: + num_static_params = 0 + num_grid_points = grid_size ** num_params + + # define the grid axes across each dimension + grid_axes = [] + for n in range(0, num_params): + step = (max_vals[n] - min_vals[n]) / (grid_size + 1) + grid_axes.append(np.linspace(offsets[n], step*(grid_size - 1) + offsets[n], grid_size)) + unconstrained_vals = c2u_array(grid_axes, min_vals, max_vals) + # transform the grid to the unconstrained space + grid = np.meshgrid(*unconstrained_vals) + concatenated_grid = [] + for grid_points in grid: + concatenated_grid.append(np.reshape(grid_points, num_grid_points)) + + grid_points = np.asarray(concatenated_grid) + + params_temp = copy.deepcopy(current_continuous) + proposed_continuous = np.empty([dim * num_params + num_static_params]) + lp_vals = np.empty([num_grid_points]) + proposals = [] + for n in range(0, num_grid_points): + grid_point = grid_points[:, n] + param_ind1 = 0 + param_ind2 = 0 + for param_num in range(0, num_params): + if num_static_params != 0: + if param_num == 0: + param_ind1 += inds[param_num] + param_ind2 += inds[param_num] + proposed_continuous[0:param_ind2 + 1] = params_temp[0:param_ind1 + 1] + else: + old_param_ind1 = copy.deepcopy(param_ind1) + old_param_ind2 = copy.deepcopy(param_ind2) + param_ind1 += inds[param_num] - inds[param_num - 1] + param_ind2 += inds[param_num] - inds[param_num - 1] + proposed_continuous[old_param_ind2:param_ind2 + 1] = params_temp[old_param_ind1:param_ind1 + 1] + proposed_continuous[param_ind2:(param_ind2 + dim - 1)] = params_temp[param_ind1:(param_ind1 + dim - 1)] + proposed_continuous[param_ind2 + dim - 1] = grid_point[param_num] + param_ind1 += dim - 1 + param_ind2 += dim # assuming that we're only adding one to the dimension + if param_ind1 == (total_params - 1): + proposed_continuous[param_ind2] = params_temp[param_ind1] + else: + proposed_continuous[param_ind2:-1] = params_temp[param_ind1:-1] + lp_vals[n] = stan_model.eval(proposed_data, proposed_continuous[0:param_length])[0] + proposals.append(proposed_continuous) + + # randomly sample grid value using the distribution of log_prob values at each point in the grid + max_lp = np.max(lp_vals) # need to normalise log_prob values in order to do sum over prob + norm_lp_vals = lp_vals - max_lp + logprobs = norm_lp_vals - logsumexp(norm_lp_vals) + + return proposals, logprobs + + +def forward_grid_search(data_function, stan_model, grid_size, min_vals, max_vals, offsets, inds, current_continuous, + proposed_discrete): + # input: + # data_function - function returning data passed to Stan model for given discrete parameters + # stan_model - initialised stan_model object + # grid_size - number of points across which to sample each new parameter + # min_vals - list of minimum values of each constrained parameter in order that they appear in + # current_continuous + # max_vals - list of maximum values of each constrained parameter in order that they appear in + # current_continuous + # offsets - randomly sampled offsets used to define the grids + # inds - indices of parameters to sample (starting with 0) + # current_continuous - current unconstrained parameters (including those with unchanging dimensions) + # proposed_discrete - proposed discrete parameter (assumed to be 1 parameter controlling birth / death moves) + + dim = proposed_discrete.value + + proposed_data = data_function(proposed_discrete) + + param_length = stan_model.num_unconstrained_parameters(proposed_data) + + num_params = len(min_vals) + + total_params = len(current_continuous) + num_static_params = total_params - (dim - 1) * num_params + if num_static_params < 0: + num_static_params = 0 + num_grid_points = grid_size ** num_params + + # define the grid axes across each dimension + grid_axes = [] + for n in range(0, num_params): + step = (max_vals[n] - min_vals[n]) / (grid_size + 1) + grid_axes.append(np.linspace(offsets[n], step*(grid_size - 1) + offsets[n], grid_size)) + unconstrained_vals = c2u_array(grid_axes, min_vals, max_vals) + # transform the grid to the unconstrained space + grid = np.meshgrid(*unconstrained_vals) + concatenated_grid = [] + for grid_points in grid: + concatenated_grid.append(np.reshape(grid_points, num_grid_points)) + + grid_points = np.asarray(concatenated_grid) + + params_temp = copy.deepcopy(current_continuous) + proposed_continuous = np.empty([dim * num_params + num_static_params]) + lp_vals = np.empty([num_grid_points]) + for n in range(0, num_grid_points): + grid_point = grid_points[:, n] + param_ind1 = 0 + param_ind2 = 0 + for param_num in range(0, num_params): + if num_static_params != 0: + if param_num == 0: + param_ind1 += inds[param_num] + param_ind2 += inds[param_num] + proposed_continuous[0:param_ind2 + 1] = params_temp[0:param_ind1 + 1] + else: + old_param_ind1 = copy.deepcopy(param_ind1) + old_param_ind2 = copy.deepcopy(param_ind2) + param_ind1 += inds[param_num] - inds[param_num - 1] + param_ind2 += inds[param_num] - inds[param_num - 1] + proposed_continuous[old_param_ind2:param_ind2 + 1] = params_temp[old_param_ind1:param_ind1 + 1] + proposed_continuous[param_ind2:(param_ind2 + dim - 1)] = params_temp[param_ind1:(param_ind1 + dim - 1)] + proposed_continuous[param_ind2 + dim - 1] = grid_point[param_num] + param_ind1 += dim - 1 + param_ind2 += dim # assuming that we're only adding one to the dimension + if param_ind1 == (total_params - 1): + proposed_continuous[param_ind2] = params_temp[param_ind1] + else: + proposed_continuous[param_ind2:-1] = params_temp[param_ind1:-1] + lp_vals[n] = stan_model.eval(proposed_data, proposed_continuous[0:param_length])[0] + + # randomly sample grid value using the distribution of log_prob values at each point in the grid + max_lp = np.max(lp_vals) # need to normalise log_prob values in order to do sum over prob + p_vals = np.exp(lp_vals - max_lp) # need to ensure this never underflows / overflows + probs = p_vals / np.sum(p_vals) + samp = np.random.choice(range(0, grid_size), p=probs) + param_ind1 = 0 + param_ind2 = 0 + for param_num in range(0, num_params): + if num_static_params != 0: + if param_num == 0: + param_ind1 += inds[param_num] + param_ind2 += inds[param_num] + proposed_continuous[0:param_ind2 + 1] = params_temp[0:param_ind1 + 1] + else: + old_param_ind1 = copy.deepcopy(param_ind1) + old_param_ind2 = copy.deepcopy(param_ind2) + param_ind1 += inds[param_num] - inds[param_num - 1] + param_ind2 += inds[param_num] - inds[param_num - 1] + proposed_continuous[old_param_ind2:param_ind2 + 1] = params_temp[old_param_ind1:param_ind1 + 1] + proposed_continuous[param_ind2:(param_ind2 + dim - 1)] = params_temp[param_ind1:(param_ind1 + dim - 1)] + proposed_continuous[param_ind2 + dim - 1] = grid_points[param_num, samp] + param_ind1 += dim - 1 + param_ind2 += dim + proposed_continuous[param_ind2:-1] = params_temp[param_ind1:-1] + + return proposed_continuous, np.log(probs[samp]) + + +def reverse_grid_search(data_function, stan_model, grid_size, min_vals, max_vals, inds, current_continuous, current_discrete, + death_ind): + # input: + # data_function - function returning data passed to Stan model for given discrete parameters + # stan_model - initialised stan_model object + # grid_size - number of points across which to sample each new parameter + # min_vals - list of minimum values of each constrained parameter in order that they appear in + # current_continuous + # max_vals - list of maximum values of each constrained parameter in order that they appear in + # current_continuous + # inds - indices of parameters to sample (starting with 0) + # current_continuous - current unconstrained parameters (including those with unchanging dimensions) + # current_discrete - current discrete parameter (assumed to be 1 parameter controlling birth / death moves) + # death_ind - index of discrete variable proposed for death + + dim = current_discrete.value + + current_data = data_function(current_discrete) + + param_length = stan_model.num_unconstrained_parameters(current_data) + + num_params = len(min_vals) + num_grid_points = grid_size ** num_params + + # define the grid axes across each dimension + grid_axes = [] + param_ind = 0 + for param_num in range(0, num_params): + if param_num == 0: + param_ind += inds[param_num] + else: + param_ind += inds[param_num] - inds[param_num - 1] + step = (max_vals[param_num] - min_vals[param_num]) / (grid_size + 1) + offset = np.mod(u2c(current_continuous[param_ind + death_ind], min_vals[param_num], max_vals[param_num]), step) + grid_axes.append(np.linspace(offset, step*(grid_size - 1) + offset, grid_size)) + param_ind += dim + + unconstrained_vals = c2u_array(grid_axes, min_vals, max_vals) + # transform the grid to the unconstrained space + grid = np.meshgrid(*unconstrained_vals) + concatenated_grid = [] + for grid_points in grid: + concatenated_grid.append(np.reshape(grid_points, num_grid_points)) + + grid_points = np.asarray(concatenated_grid) + + proposed_continuous = copy.deepcopy(current_continuous) + lp_vals = np.empty([num_grid_points]) + for n in range(0, num_grid_points): + grid_point = grid_points[:, n] + param_ind1 = 0 + for param_num in range(0, num_params): + if param_num == 0: + param_ind1 += inds[param_num] + else: + param_ind1 += inds[param_num] - inds[param_num - 1] + param_ind2 = param_ind1 + death_ind + proposed_continuous[param_ind2] = grid_point[param_num] + param_ind1 += dim - 1 + lp_vals[n] = stan_model.eval(current_data, proposed_continuous[0:param_length])[0] + + # calculate probability of choosing parameter values from specified grid + max_lp = np.max(lp_vals) # need to normalise log_prob values in order to do sum over prob + p_vals = np.exp(lp_vals - max_lp) # need to ensure this never underflows / overflows + probs = p_vals / np.sum(p_vals) + grid_point = np.empty([num_params]) + for param_num in range(0, num_params): + if param_num == 0: + param_ind += inds[param_num] + else: + param_ind += inds[param_num] - inds[param_num - 1] + grid_point[param_num] = current_continuous[param_ind] + param_ind += dim - 1 + ind = (np.abs(grid_points - grid_point)).argmin() + return np.log(probs[ind]) diff --git a/discretesampling/base/algorithms/rjmcmc.py b/discretesampling/base/algorithms/rjmcmc.py new file mode 100644 index 0000000..f5a2ea9 --- /dev/null +++ b/discretesampling/base/algorithms/rjmcmc.py @@ -0,0 +1,154 @@ +import math +import copy +import numpy as np +from scipy.special import logsumexp +from discretesampling.base.random import RNG +from discretesampling.base.algorithms.continuous.NUTS import NUTS +from discretesampling.base.algorithms.continuous.RandomWalk import RandomWalk + + +class DiscreteVariableRJMCMC(): + + def __init__(self, variableType, discrete_target, + stan_model, data_function, continuous_proposal_type, continuous_update, accept_reject=False, + always_update=False, do_warmup=True, update_probability=0.5, initialProposal=None, warmup_iters=100): + + self.variableType = variableType + self.proposalType = variableType.getProposalType() + self.initialProposal = initialProposal + self.discrete_target = discrete_target + + self.data_function = data_function + + self.data_function = data_function + self.continuous_proposal_type = continuous_proposal_type + self.continuous_update = continuous_update + self.stan_model = stan_model + self.accept_reject = accept_reject + self.always_update = always_update + self.update_probability = update_probability + + self.rng = RNG() + + # initialise samplers for continuous parameters + if self.continuous_update is RandomWalk: + self.csampler = RandomWalk(self.stan_model, self.data_function, self.rng) + elif self.continuous_update is NUTS: + if do_warmup: + self.csampler = NUTS(True, self.stan_model, self.data_function, self.rng, 0.9, warmup_iters) + else: + self.csampler = NUTS(False, self.stan_model, self.data_function, self.rng, 0.9, warmup_iters) + else: + raise NameError("Continuous update type not defined") + + def propose(self, current_discrete, current_continuous): + proposed_discrete = current_discrete + q = self.rng.random() + r0 = 0 + r1 = 0 + # Update vs Birth/death + if (q < self.update_probability): + # Perform update in continuous space + [proposed_continuous, r0, r1] = self.csampler.sample(current_continuous, current_discrete) + else: + # Perform discrete update + forward_proposal = self.proposalType(current_discrete) + proposed_discrete = forward_proposal.sample() + + # Birth/death continuous dimensions + # It would probably be better if this was called as above, with separate eval() and sample() functions. + # However, the random choices made in the forward move would need to be passed to calculate the probability of + # the reverse move (e.g. random selection for birth / death). This information might be better stored in the + # discrete variables. + if self.always_update: + # need to do this to initialise NUTS (if we don't do this then we can end up in a situation in SMC where we're + # comparing NUTS proposals between two starting points without sampled momenta / NUTS parameters) + init_proposed_continuous = self.continuous_proposal_type().sample(current_discrete, current_continuous, + proposed_discrete, + self.rng) + [proposed_continuous, r0, r1] = self.csampler.sample(init_proposed_continuous, proposed_discrete) + else: + proposed_continuous = self.continuous_proposal.sample(current_discrete, current_continuous, proposed_discrete, + self.rng) + + return proposed_discrete, proposed_continuous, r0, r1 + + def sample(self, N): + current_discrete = self.initialProposal.sample() + self.continuous_proposal = self.continuous_proposal_type() + if hasattr(current_discrete.value, "__len__"): + empty_discrete = self.variableType(np.zeros(current_discrete.value.shape())) + else: + empty_discrete = self.variableType(0) + + # get initial continuous proposal by performing a birth move from a 0-dimensional model + current_continuous = self.continuous_proposal.sample(empty_discrete, np.array([]), current_discrete, self.rng) + samples = [] + for i in range(N): + [proposed_discrete, proposed_continuous, r0, r1] = self.propose(current_discrete, current_continuous) + log_acceptance_ratio = self.eval(current_discrete, current_continuous, proposed_discrete, proposed_continuous, + r0, r1) + if log_acceptance_ratio > 0: + log_acceptance_ratio = 0 + acceptance_probability = min(1, math.exp(log_acceptance_ratio)) + q = self.rng.random() + # Accept/Reject + if (q < acceptance_probability): + current_discrete = proposed_discrete + current_continuous = proposed_continuous + print("Iteration {}, params = {}".format(i, current_continuous)) + samples.append([copy.copy(current_discrete), current_continuous]) + + return samples + + def eval(self, current_discrete, current_continuous, proposed_discrete, proposed_continuous, r0, r1): + if current_discrete == proposed_discrete: + # continuous parameter update move + if self.continuous_update == NUTS and self.accept_reject is False: + log_acceptance_ratio = 0 + else: + log_acceptance_ratio = self.csampler.eval(current_continuous, current_discrete, proposed_continuous, r0, r1) \ + - np.log(self.update_probability) + else: + # discrete move + forward_proposal = self.proposalType(current_discrete) + reverse_proposal = self.proposalType(proposed_discrete) + discrete_forward_logprob = forward_proposal.eval(proposed_discrete) + discrete_reverse_logprob = reverse_proposal.eval(current_discrete) + + if self.always_update and not (self.continuous_update == NUTS and self.accept_reject is False): + # we need to sum over all of the possible combinations of birth / death moves + updates to get to + # proposed_continuous + continuous_proposals = self.continuous_proposal_type().eval_all(current_discrete, current_continuous, + proposed_discrete, proposed_continuous) + continuous_proposal_logprobs = [] + for continuous_proposal, continuous_proposal_logprob in continuous_proposals: + nuts_logprob = self.csampler.eval(continuous_proposal, proposed_discrete, proposed_continuous, r0, r1) + continuous_proposal_logprobs.append(continuous_proposal_logprob + nuts_logprob) + continuous_proposal_logprob = logsumexp(continuous_proposal_logprobs) + else: + continuous_proposal_logprob = self.continuous_proposal.eval(current_discrete, current_continuous, + proposed_discrete, proposed_continuous) + + # Setup data for stan model + current_data = self.data_function(current_discrete) + proposed_data = self.data_function(proposed_discrete) + param_length = self.stan_model.num_unconstrained_parameters(current_data) + p_param_length = self.stan_model.num_unconstrained_parameters(proposed_data) + + # P(theta | discrete_variables) + current_continuous_target = self.stan_model.eval(current_data, current_continuous[0:param_length])[0] + proposed_continuous_target = self.stan_model.eval(proposed_data, proposed_continuous[0:p_param_length])[0] + + # P(discrete_variables) + current_discrete_target = self.discrete_target.eval(current_discrete) + proposed_discrete_target = self.discrete_target.eval(proposed_discrete) + + jacobian = 0 # math.log(1) + + log_acceptance_ratio = (proposed_continuous_target - current_continuous_target + + proposed_discrete_target - current_discrete_target + + discrete_reverse_logprob - discrete_forward_logprob + + continuous_proposal_logprob + jacobian - np.log(1 - self.update_probability)) + + return log_acceptance_ratio diff --git a/discretesampling/base/random.py b/discretesampling/base/random.py index 2482760..66242ed 100644 --- a/discretesampling/base/random.py +++ b/discretesampling/base/random.py @@ -25,3 +25,9 @@ def randomChoice(self, choices): def randomChoices(self, population, weights=None, cum_weights=None, k=1): return self.nprng.choice(population, size=k, replace=True, p=weights) + + def randomMvNormal(self, Mu, Sigma): + return self.nprng.multivariate_normal(Mu, Sigma) + + def randomNormal(self, mu=0, sigma=1): + return self.nprng.normal(mu, sigma) diff --git a/discretesampling/base/reversible_jump.py b/discretesampling/base/reversible_jump.py new file mode 100644 index 0000000..8ea529b --- /dev/null +++ b/discretesampling/base/reversible_jump.py @@ -0,0 +1,170 @@ +from ..base import types +from ..base.algorithms import rjmcmc +from discretesampling.base.random import RNG +from discretesampling.base.algorithms.continuous import NUTS +import numpy as np +import copy + + +def set_proposal_attributes( + base_DiscreteVariableTarget, + stan_model, + data_function, + continuous_proposal_type, + continuous_sampler_type, + update_probability=0.5 +): + global base_target + global model + global data_func + global cont_proposal_type + global cont_sampler_type + global update_prob + + # attributes used in proposal (must be set prior to sampling) + base_target = base_DiscreteVariableTarget + model = stan_model + data_func = data_function + cont_proposal_type = continuous_proposal_type + cont_sampler_type = continuous_sampler_type + update_prob = update_probability + + +# ReversibleJumpVariable inherits from DiscreteVariable +class ReversibleJumpVariable(types.DiscreteVariable): + def __init__(self, DiscreteVariables, ContinuousVariables): + self.discrete = DiscreteVariables + self.continuous = ContinuousVariables + + # need to store the current and proposed momenta, plus NUTS parameters if we are using NUTS + self.r0 = 0 + self.r1 = 0 + self.NUTS_params = {} + + @classmethod + def getProposalType(self): + return ReversibleJumpProposal + + @classmethod + def getTargetType(self): + return ReversibleJumpTarget + + +# ReversibleJumpProposal uses discrete proposals defined by other DiscreteVariableProposal +# The main purpose of this proposal class proposes new continuous parameters based on the discrete proposal +class ReversibleJumpProposal(types.DiscreteVariableProposal): + def __init__(self, currentReversibleJumpVariable, rng=RNG()): + self.currentReversibleJumpVariable = currentReversibleJumpVariable + self.rjmcmc = rjmcmc.DiscreteVariableRJMCMC( + type(currentReversibleJumpVariable.discrete), + base_target, + model, + data_func, + cont_proposal_type, + cont_sampler_type, + True, True, False, + update_prob + ) + # copy over previous NUTS parameters + if cont_sampler_type == NUTS: + self.rjmcmc.csampler.NUTS_params = self.currentReversibleJumpVariable.NUTS_params + + self.rng = rng + + @classmethod + def norm(self, x): + return self.base_DiscreteVariableProposal.norm(x.discrete) + + @classmethod + def heuristic(self, x, y): + # Proposal can be type specified by discrete variables or continuous update + return self.base_DiscreteVariableProposal.heuristic(x.discrete, y.discrete) or x.discrete == y.discrete + + def sample(self, target=None): + + proposedReversibleJumpVariable = copy.deepcopy(self.currentReversibleJumpVariable) + + current_discrete = self.currentReversibleJumpVariable.discrete + current_continuous = self.currentReversibleJumpVariable.continuous + + [proposed_discrete, proposed_continuous, r0, r1] = self.rjmcmc.propose(current_discrete, current_continuous) + proposedReversibleJumpVariable.discrete = proposed_discrete + proposedReversibleJumpVariable.continuous = proposed_continuous + + if cont_sampler_type == NUTS: + proposedReversibleJumpVariable.r0 = r0 + proposedReversibleJumpVariable.r1 = r1 + proposedReversibleJumpVariable.NUTS_params = self.rjmcmc.csampler.NUTS_params + + return proposedReversibleJumpVariable + + def eval(self, proposedReversibleJumpVariable): + current_discrete = self.currentReversibleJumpVariable.discrete + current_continuous = self.currentReversibleJumpVariable.continuous + + proposed_discrete = proposedReversibleJumpVariable.discrete + proposed_continuous = proposedReversibleJumpVariable.continuous + + r0 = proposedReversibleJumpVariable.r0 + r1 = proposedReversibleJumpVariable.r1 + + log_acceptance_ratio = self.rjmcmc.eval(current_discrete, current_continuous, proposed_discrete, proposed_continuous, + r0, r1) + + if log_acceptance_ratio > 0: + log_acceptance_ratio = 0 + + return log_acceptance_ratio + + +class ReversibleJumpTarget(types.DiscreteVariableTarget): + + def eval(self, proposedReversibleJumpVariable): + proposed_discrete = proposedReversibleJumpVariable.discrete + proposed_continuous = proposedReversibleJumpVariable.continuous + proposed_data = data_func(proposed_discrete) + param_length = model.num_unconstrained_parameters(proposed_data) + proposed_continuous_target = model.eval(proposed_data, proposed_continuous[0:param_length])[0] + proposed_discrete_target = base_target.eval(proposed_discrete) + target = proposed_continuous_target + proposed_discrete_target + return target + + +class ReversibleJumpInitialProposal(types.DiscreteVariableInitialProposal): + def __init__(self, base_DiscreteVariableType, base_DiscreteVariableInitialProposal): + self.base_type = base_DiscreteVariableType + self.base_proposal = base_DiscreteVariableInitialProposal + + def sample(self, rng=RNG(), target=None): + rjmcmc_proposal = rjmcmc.DiscreteVariableRJMCMC(self.base_type, base_target, model, + data_func, cont_proposal_type, cont_sampler_type, True, True, False, update_prob) + proposed_discrete = self.base_proposal.sample() + if hasattr(proposed_discrete.value, "__len__"): + empty_discrete = self.base_type(np.zeros(proposed_discrete.value.shape())) + else: + empty_discrete = self.base_type(0) + # get initial continuous proposal by performing a birth move from a 0-dimensional model + proposed_continuous = cont_proposal_type().sample(empty_discrete, np.array([]), proposed_discrete, rng=rng) + + # do an initial continuous move in case we need to initialise NUTS + [proposed_continuous, r0, r1] = rjmcmc_proposal.csampler.sample(proposed_continuous, proposed_discrete) + + proposedReversibleJumpVariable = ReversibleJumpVariable(proposed_discrete, proposed_continuous) + + if cont_sampler_type == NUTS: + proposedReversibleJumpVariable.r0 = r0 + proposedReversibleJumpVariable.r1 = r1 + proposedReversibleJumpVariable.NUTS_params = rjmcmc_proposal.csampler.NUTS_params + + return proposedReversibleJumpVariable + + def eval(self, proposedReversibleJumpVariable, target=None): + proposed_discrete = proposedReversibleJumpVariable.discrete + if hasattr(proposed_discrete.value, "__len__"): + empty_discrete = self.base_type(np.zeros(proposed_discrete.value.shape())) + else: + empty_discrete = self.base_type(0) + proposed_continuous = proposedReversibleJumpVariable.continuous + continuous_logprob = cont_proposal_type().eval(empty_discrete, np.array([]), proposed_discrete, proposed_continuous) + discrete_logprob = self.base_proposal.eval(proposed_discrete) + return continuous_logprob + discrete_logprob diff --git a/discretesampling/base/stan_model.py b/discretesampling/base/stan_model.py new file mode 100644 index 0000000..cb23331 --- /dev/null +++ b/discretesampling/base/stan_model.py @@ -0,0 +1,60 @@ +import math +import os +import logging +import bridgestan as bs +from cmdstanpy import write_stan_json +import numpy as np + + +class StanModel: + def __init__(self, model_file): + self.model_file = os.path.abspath(model_file) + self.model_filename = os.path.basename(self.model_file) + self.model_path = os.path.dirname(self.model_file) + self.exec_name = self.model_filename.replace(".stan", "") + self.exec_path = os.path.join(self.model_path, self.exec_name) + self.data_file = self.exec_name + ".data.json" + self.compiled = False + self.model = None + self.lib = None + self.data = None + + def compile(self): + if not self.compiled: + logging.info("Compiling Stan model ", self.exec_name, "with bridge-stan...") + self.model = bs.StanModel.from_stan_file(self.model_file, self.data_file) + logging.info("Finished compiling Stan model ", self.exec_name, ".") + self.compiled = True + + def num_unconstrained_parameters(self, data): + self.prepare_data(data) + self.compile() + return self.model.param_unc_num() + + def eval(self, data, params): + self.prepare_data(data) + self.compile() + + logprob = -math.inf + gradients = None + + try: + if self.model.param_unc_num() != len(params): + raise ValueError("Array of incorrect length passed to log_density;" + + " expected {x}, found {y}".format(x=self.model.param_unc_num(), y=len(params))) + + logprob, gradients = self.model.log_density_gradient(theta_unc=np.array(params)) + + except RuntimeError: + logging.error("Something went wrong when trying to evaluate the stan model") + + return logprob, gradients + + def prepare_data(self, data): + if self.data != data: + # Cache data + self.data = data + # Write params to data file + write_stan_json(self.data_file, dict(data)) + # If the data changes, the model also needs to be recompiled before it's used + self.compiled = False diff --git a/discretesampling/base/util.py b/discretesampling/base/util.py index af96e0e..0a948b6 100644 --- a/discretesampling/base/util.py +++ b/discretesampling/base/util.py @@ -1,4 +1,5 @@ import numpy as np +from scipy.special import logit, expit def pad(x, exec): @@ -48,3 +49,29 @@ def gather_all(particles, exec): all_particles = restore(all_x, all_particles) return all_particles + +# transforms a continuous unconstrained parameter to the constrained space for a given set of limit values +def u2c(unconstrained_param, min_val, max_val): + return min_val+(max_val-min_val)*expit(unconstrained_param) + + +# transforms a continuous constrained parameter to the unconstrained space for a given set of limit values +def c2u(constrained_param, min_val, max_val): + return logit((constrained_param - min_val) / (max_val - min_val)) + + +# transforms arrays of continuous unconstrained parameters to constrained parameters for a given set of limit values +def u2c_array(unconstrained_params, min_vals, max_vals): + constrained_params = [] + for param_array, min, max in zip(unconstrained_params, min_vals, max_vals): + constrained_params.append(min+(max-min)*expit(param_array)) + return constrained_params + + +# transforms arrays of continuous constrained parameters to unconstrained parameters for a given set of limit values +def c2u_array(constrained_params, min_vals, max_vals): + unconstrained_params = [] + for param_array, min, max in zip(constrained_params, min_vals, max_vals): + unconstrained_params.append(logit((param_array-min)/(max-min))) + return unconstrained_params + diff --git a/discretesampling/domain/additive_structure/additive_target.py b/discretesampling/domain/additive_structure/additive_target.py index ffbd903..24ebb3e 100644 --- a/discretesampling/domain/additive_structure/additive_target.py +++ b/discretesampling/domain/additive_structure/additive_target.py @@ -16,11 +16,13 @@ def eval(self, x): # presumably some function of x.discrete_set and some data which # could be defined in constructor as self.data - y = [self.evaluate(x.discrete_set, self.xtrain.iloc[i]) for i in range(self.xtrain.shape[0])] + y = [self.evaluate(x.discrete_set, self.xtrain.iloc[i]) + for i in range(self.xtrain.shape[0])] cov = 0.25 # calculate Π(Y_i|M,theta,x_i) - target1 = [self.log_likelihood(self.ytrain.iloc[i], y[i], cov) for i in range(len(y))] + target1 = [self.log_likelihood( + self.ytrain.iloc[i], y[i], cov) for i in range(len(y))] targ = sum(target1) # p(M) diff --git a/examples/5_targets_noisy.data.json b/examples/5_targets_noisy.data.json new file mode 100644 index 0000000..cf9fee2 --- /dev/null +++ b/examples/5_targets_noisy.data.json @@ -0,0 +1,2 @@ +{"L":25,"T":200,"fs":2000,"wave_speed":1481,"d":[[0.61534684304992771,0.83642894222758457,1.8181178278536887,2.2615886858554362,3.3324092987564589,1.3211999258237543,2.8918868502650996,3.2642896832882675,3.04114951536648,3.7508072575449574,4.9169358431489592,4.3639314565484533,4.8766068587368947,5.2101968238686229,4.1805790558842428,3.99379855221554,3.3573920634556416,3.6446524180764284,3.1497440825431307,3.0848362233596456,2.7403779904623633,2.1270414808605929,1.5445252285496709,1.5494511944832614,1.3707956611439942,0.55216436259507806,0.83096040647117664,0.26220440490748465,-0.086646630882627251,0.24336890618767104,0.72218171512341489,-0.18549199369463198,-0.42527407825572083,-0.84052482785478821,-0.83439384622683432,-0.62735086187346911,-1.1789271776279464,-1.6027676101259241,-0.57048953739005093,-1.0444608968256022,-0.39092675629024992,-0.88919913804339057,-0.9848581786595999,-0.94826479232032379,-0.0044100627925824087,0.91848710860332894,0.10738992732671378,0.29405875087895461,0.13477628181825058,0.16209121436290208,1.1827375446153692,-0.46917187475541205,-0.2672393390880039,-0.66281031376269972,-0.1331847732863139,-0.54922186198475276,-1.0448410427688208,-0.72705476993052676,-1.8245089078588719,-1.4858712569097283,-1.663134878733751,-1.1814501320280812,-1.8848258890776788,-2.2978951138899597,-3.1143365618514753,-3.2715536902958697,-2.2926274776237503,-2.9455900198512155,-2.6768689735356332,-2.9035751684883517,-2.2422209369328421,-2.0357241388443121,-1.0485329683624294,-0.74102392896648162,-0.81039628538256014,-0.45800730720150773,-1.2174067520293743,0.39195218217759542,-0.31849407838026905,-0.24566804765541894,2.5341830726984478,1.4942596241035635,1.3659680009872945,2.0346503874468036,2.36180074548961,1.5455592328958077,2.4269120362049685,3.1994656855474664,3.0449901871693252,2.2279893172258065,2.7845025715159504,2.8647385861252039,2.8301197436483,1.7426593784730853,2.4492310695841852,2.7777032392674168,0.86676642165818008,1.4979770418820957,1.5392890928864706,0.62846889133009709,0.29453212741734797,-0.67581414127880879,-0.96005790170145122,-1.8509396087234504,-1.7238827535208681,-2.2894865553945039,-2.6505794138357586,-1.8584223651057892,-2.847982252659591,-2.6248046304599124,-2.8425951384061583,-3.5668804132958338,-3.9721491395867594,-3.6516842048925549,-2.6467165568736704,-2.7311772618705441,-1.9475445814973713,-2.2206611826084175,-2.358587804153065,-1.2867601795175494,-0.40053255833957857,-0.421347744080685,-0.24526980023545594,0.03546797013108885,1.2666318046405631,1.4571535782238707,1.1263870099252196,0.97676407241296836,1.4641177532244809,2.3163274081279455,1.7867281443398548,2.2753911301231118,2.495056029087428,2.2656094489563863,2.2757042337587272,2.4229432455496425,2.1680430448642722,1.6432929261836562,2.5224343020161064,1.508601314320269,2.2808754043670256,1.3738779103421632,1.1677155113731066,1.7458937296384685,1.0117541777446559,0.475680279619593,0.0767310467889441,0.0673470989798759,0.38774938590457253,0.50210488410412957,0.024005478351137689,-0.29350480925127942,0.1048467728925758,0.28680023271700367,-0.21800208028324422,0.78034574300978021,-0.03717970190179401,0.69569380902695066,0.84232150605655032,0.5022417187680297,0.54545201933363141,1.1295329773445935,-0.11678219767893561,0.87335825096642172,0.93700634089551738,1.1008588873821228,0.45685725483090245,0.13779716716381496,-0.57257028980491187,0.88186455418091858,0.96991266655503616,-0.0096674418165365039,1.1459819226397183,-0.65475611782235732,-1.3104815138598014,-1.0728440305484519,-1.3105901914460187,-1.3971187015231843,-2.0563467228843395,-2.8730340559869743,-3.6146810847317963,-3.48782540585491,-4.188962416215845,-3.421503404532257,-4.0677656755102634,-5.1251101470544356,-4.1429162163800513,-4.8539472671259816,-3.7211254766243105,-4.3081857916106134,-4.0912026066256537,-4.4051152837587653,-3.4914527170305472,-3.4011549278310276,-2.9497285653605649,-2.7913482525237194,-1.8388194957381276,-2.8729475969411293,-1.0352896775571208,-0.039772798670180354],[-0.47590026639917238,0.21874527930489329,-0.077311328059603224,0.57402425425967385,0.43596220557165877,1.511228402090302,1.3530817026220723,1.953073540465724,1.3586207320409809,1.4322045775919969,1.6701687905325593,2.2606971755661065,0.33839351450423871,0.66932673348240768,2.0959146129821096,2.3489927874856593,1.2284458887274141,1.9295244303618957,2.2569051320178293,1.187333179732029,0.79070296335384471,0.714607765672214,0.77162743797034,1.2597168523089852,0.75178902980178508,1.8497347696770059,0.67294675655101388,0.51411422043970312,1.1902239164657515,1.0966654800965721,1.2842828945464395,0.090339349483503639,1.1535880664374205,0.40661087082183356,1.5104323865531613,1.8902090323299199,2.1789146551573637,2.0344342502415458,1.9523606150057911,2.1303804674079312,1.7234595173699976,1.71577352011669,2.2004874610349363,2.1597056237731218,1.2762750736941308,0.47851875228660468,0.065591614644541663,0.21023628113520121,0.30424143433598805,0.78807753192775065,-0.082008711933017453,-0.0099718734173353737,-0.0657095147318662,-1.2313862890650809,-1.250680554534839,-1.4813980380270595,-1.7770905814871396,-1.8482190892654478,-1.4651271949487485,-2.3154793978236019,-1.8861838874245918,-2.1995939502423649,-2.372010767805889,-2.0380707589378755,-1.4711520452635429,-1.7806081187457568,-3.2370214443001686,-2.0696524398341789,-2.8085701036927357,-3.0049203505302415,-2.3867630021356683,-2.6261420714176849,-1.831478845163673,-0.30171303528235582,-0.91835684783494032,0.96464622536683164,0.12350913161636659,0.72711999840692243,0.43629936838070216,0.63507952451444,0.62481961601855729,1.3098142839398697,0.66795826613294351,1.128115732588044,1.1233572866150983,1.8063147405957942,2.5027836768533995,2.8048630212418209,0.789063711784665,1.7425029675396819,1.1083671593751714,0.72847795639380331,1.1263294461511708,0.62044200976364261,0.64513474773105606,0.7539913306491427,-0.355168390395348,-0.22249195902609023,-1.2508327435749704,-0.83705994967159325,-0.38081464261281794,-1.9581185210110152,-1.0758244226519902,-2.1055394504315452,-1.4477013139405304,-1.8771948990458469,-1.85197883326378,-2.5405827427282492,-2.1856183549786188,-1.7846248367578197,-1.5787109732584583,-1.6075721009588935,-1.7039728010587192,-0.56828630275232683,-1.160884498530325,-1.2690177810734962,-1.0687860053752933,-0.5313637067479563,0.16643428169468516,-0.31978917034151394,0.50964791983379232,0.35542995942144329,0.11234059487902409,0.11792246492173264,0.72621704594113778,0.67383829842737841,0.337080600632408,0.43248876101703204,1.1159255387936655,1.1018076838079209,0.65974336099462172,1.1853532433604914,0.539574676967853,0.58790813135898867,1.1868088140073381,0.21977142794734145,-0.61075397877559556,-0.30803161964089709,0.41610050409659616,0.10926122458009092,1.554785521244892,-0.54819085363212872,0.03816336333385438,0.20311064633917786,-0.45796857042889572,0.83778922386585564,0.76217281224399192,1.2751428255744912,0.98048021397027751,-0.988534872526046,0.4663691682477043,1.2787907370194243,0.893684142514213,1.6218553897029466,1.6216844947843898,1.5462924237832403,2.0048282132190374,2.3586094506843658,2.7808483459054303,1.8410262979979142,2.0484865701542336,1.9848495271315603,2.741777246056861,3.0368603194630124,2.6702077513178222,2.6636252324111829,1.9282718599309427,1.1747021084106062,1.6512811764283728,1.8777323350495969,0.94210671310522676,0.80441921307242337,-0.79960338460413682,-0.12252403569032011,-0.87150529031992985,-0.9042029752756342,-0.92207171981091551,-2.2465298514584719,-1.5206544282333976,-3.0665512315808154,-2.8645510915424621,-2.9854951167648571,-2.8191487309167158,-3.6837552275800984,-4.1130676094179961,-3.0112849414573533,-3.5633601743124443,-4.6772963529202,-3.4879515511038908,-4.3792417224149141,-2.970865456581965,-2.4700054002800163,-3.0163898593600389,-3.1096198585174819,-1.7467467465485371,-2.3004442856899474,-2.1252680450975583,-2.0099866811034981,-1.7741664155362398,-0.56774476259098439],[-1.427405877305286,-0.842859329189527,-1.23207865195116,-1.2665723774586028,-1.9722549581616502,-1.4413858099388426,-0.936557151425962,-1.1316897025040307,-0.20348527768350155,-1.6755443527460172,-0.76560596254051561,-0.53021180044924621,-0.85945322119573508,-0.352289789134884,0.51214975481137848,-0.75493859917060357,-0.31834554496272993,-0.76826619168716581,0.38608614314931644,-0.50398172187980794,-0.25751453999910134,0.9508376691401601,0.71383538748939934,-0.018253776591416415,0.92041615144439182,0.12988662794689626,0.554781558737701,1.4565812305430894,1.701235438346818,0.3803514718033586,0.61265311024162583,1.169865254069137,0.97415682496517109,0.57188121718651652,0.60003462742609015,0.32796531449384408,-0.40156486656085255,0.6701284048686984,0.16341281878122726,0.66233956206790323,2.0246161431061536,0.82179864527332069,1.5741924392284377,0.27869295890332324,1.4862885745384506,0.025682157513093329,0.735683351464385,1.032926499317437,1.2077788950313459,0.843542749462143,0.36975119111364141,0.30392631176621565,0.54977933460039519,0.39506700893312052,0.60414492949704912,0.40231263799048717,1.4770093535512145,0.96361534735775822,0.93717907478231965,0.33679035065030366,0.849613871293782,0.82187413965077938,0.1359887147437987,-0.052088681358730238,-0.49119183336753075,0.84228782606612906,-0.17410313896240226,0.27858217548196856,0.2179288094700306,0.14185808362603852,-0.30343330215346259,-0.17770842264977682,-0.29698138516018718,0.45428892351427852,-0.87277227697099236,-0.46831639805487613,0.10131244222677283,-0.28316492461278286,-1.0944206683860231,-0.20478883517512492,-1.9974936199030173,-0.93208467310559118,0.064914982024518619,-0.49691246841475523,-0.371704565007141,-0.82979763656991024,-0.98404289487307173,-0.45182496933107508,0.35716345286246431,-0.93343579757340756,-0.075470916841136271,-1.0558617477306429,-0.99157275054422589,-0.44019334536548371,-1.2048814053866193,-0.27080997877363105,-1.5719435564281019,-1.9297181449208465,-0.78883046135784218,-0.53617986899170833,-0.055532845388939656,-0.62741726581823154,-0.20539778914663648,-0.55888128512315072,-0.33690148864817676,-1.1591274764788639,-0.71123090638189079,-0.96992187258819862,-0.78670509317242432,0.48431620157978494,-1.3599087245392236,-0.63198905170637021,-0.79761245512422718,-0.43183227473148394,0.29469453038246152,-0.12091012678553958,-0.87586238984242026,-1.2292418736103183,-0.04772588488412266,-0.92145570780908814,-0.73102336965982673,-0.083156990741570208,-1.0141802831613669,0.59300753478196477,-0.061361868685350963,-0.15353178207595597,0.44699829535717972,0.52831851621206738,0.084503624049632353,0.84502983660680742,0.56230040900302258,1.0125259985248154,0.67633034046089469,0.52840961007181864,1.5741265309211068,0.3795725933110165,1.1540192301107248,1.9091863413440089,0.7454703577236752,1.3167792192792172,1.199352351169821,0.89811282447834828,1.0955265037246191,1.4236508864654303,2.2889316188838373,1.5705169119180193,1.0008483059444484,0.12528474213573981,1.2191653232359343,1.202281992181566,0.64816392525993316,1.5162551713672152,1.170922379393287,1.2296906628559672,0.65139059530330246,0.50889456402517441,0.57210629308747463,0.20469629592417937,0.45250318545098306,1.0309768776580988,-0.50604378017850371,-0.30977704663755268,0.54471819981486025,0.65156314993302322,-0.030561938323557525,0.66914155373937034,0.98785022853082116,0.65774275609684607,-0.50846560689519649,-0.13518061842999179,-0.24076499647914545,-0.02837269868600853,-0.24128740508215657,-0.37492594826647041,0.950265086035284,-0.22407182574694468,-0.18917841488979598,-0.027729721420308245,-1.0975378689632718,-0.865242688451926,-0.12479532397961238,-0.25135993760835912,-1.0765444667714685,-0.0033762735779434561,-0.67038394885729979,-0.63698874216446044,-0.72332577167154233,-1.1216489053659753,-1.8287657625112383,-1.0419732601225062,-1.8185465802232563,-1.173214171376,-1.0365221490939094,-1.0842665719780464,-1.9124749481620906,-1.0645444023787931,-1.8734602677892331,-1.1542153232187358,-0.98139221163159029,-0.47764340585367815],[-1.7583952408845596,-0.16432926506457662,-0.49095118068533822,-0.26836679705729682,-0.32799954677707693,0.25341693774340668,0.18379718548760046,0.16861136207233435,0.3564145420138411,1.0595068799814469,1.0963395421473217,1.4697167779491203,1.1364176740099263,0.88194820285564424,0.34871003186583882,1.0009119496391927,0.57900602447487937,0.69647953883155267,0.84350694950738825,-0.049321146240426439,-0.27861082007587162,-0.42225390850228167,-0.45662745945259164,0.17618683330687812,-1.5608033759786777,-0.73259701393375176,-1.7779588350744004,-1.4718934991220149,-2.1181904221125487,-2.6130393787218651,-2.63873516143447,-2.1535080103255946,-2.0800521230136955,-2.2670908949198929,-3.0152248607297687,-3.3140306816119445,-2.3559875032729947,-2.1308152009821395,-3.4431398309524384,-1.9235765472023612,-0.5975364413487485,-1.4148338693599181,-1.2893066296393199,-1.0952265039491871,-1.1052185594776396,-0.044177131375442047,0.75297648140522944,0.50620640637184411,1.5078510965500806,2.1382977992119447,1.7847210604494192,2.6325909615479057,2.9373268614920347,2.8599757769553737,4.0417181150333823,2.5624170497720495,3.4225797068690644,2.8593288857380252,2.9512424278340865,3.4754345105062718,3.7001936687423345,3.0970089527881175,3.4862140703806643,2.159044331925855,2.5574011466166535,1.8667776687653022,1.4831899619987707,0.26762626024169478,1.6830344528789034,1.9638148866770668,1.0091401000266895,1.1950036363953767,0.30297226170030955,-0.55212895918824645,-0.5575230277817893,-0.67870961727526358,0.14737632989442645,-0.74028424913031954,-0.52092352362619887,-1.271255443299435,-0.707794551318686,-1.4171024309877245,-0.92303822955481751,-0.14345003290418762,0.39216885502850152,-0.15215918174112814,-0.27606228716526054,0.27973681035027625,0.16779093637042952,-1.3715284066972435,-1.0093416121535763,-0.51582076054661563,-0.54356101334145324,-0.93871767570491549,-0.53523173957411418,-0.74545103835371329,0.16049103638070394,-0.82924757283821493,-1.3198336103288888,-1.0676901180883192,-1.3964387546926256,-1.1297817269834947,-1.59120481931194,-2.4356545669297787,-1.6172294643130132,-2.2953721505583964,-2.6052441669601332,-2.9032177215405031,-2.2274133920232781,-2.7958964963292172,-3.6005700192241461,-1.7834299020767974,-3.04330298716776,-2.6973976150012624,-2.5716236251763291,-3.2738075626250396,-1.8376222346264393,-2.1081564997800926,-1.9079781902007444,-1.0147547032635984,-1.2894361684330164,-1.8828911101485089,-0.66226433012292529,0.31066777843064264,0.72201196558275538,0.63666900274703375,0.039252306826008776,1.8018800554228456,1.5718729330433368,2.0162658633156267,2.8477018629914657,2.6569903303979605,0.80044251561636259,1.9894778922668039,2.518434415667262,2.3476311045058216,3.0800993625526689,3.0392542759028527,2.5280566603774108,3.324735054512244,2.1689806714330597,2.1752939386250287,1.7667108948961272,1.3775431306228447,1.5578493152541901,0.52429964050353439,0.74906944826165234,1.0039871991184177,0.68605165000711466,-0.99033025610416037,-0.091061238386811461,-0.77565873817010511,-0.54388082100599888,-0.87171144066859729,-1.26734545379253,-1.2281194087058507,-0.98897276585363536,-0.96036202057635611,-1.688327976132443,-1.9285585713055586,-1.3896350864921685,-0.5525385086314093,-0.74487339708962585,-0.5490790727780881,-0.77844734606327359,-0.73802983086996743,0.89254281937045921,0.012004567622936148,0.25402072782625562,0.32973395614312895,0.682001984519869,0.64780893184543309,0.42156948815939144,1.2985011774636228,1.8498419052223336,1.2680405291024264,0.81599080962743908,1.195500284121908,1.624275043743675,1.9218667797340956,0.21570983449343761,0.83850712133479621,0.40708124585233957,0.72175257469543419,0.3480216159458896,-0.22309323948220339,-1.0258586430459773,-0.91866300524928346,-0.79473020345657164,-0.79407228277518749,-0.21627555306825919,-1.2256179242544325,-0.8404155659398298,-0.64291741370394062,-1.01886569536862,-1.0996587836741389,-1.4230813193505059,-0.64206530410652563,-0.16414543561099326,-0.62969575077767859],[1.5408383276030189,1.83249987434605,2.5366104530908231,2.1283723816591511,2.2742478382728026,2.9886354929531067,2.7120634586697103,3.3994660103477958,2.2520488219838146,3.6175175607493428,2.4983205259268266,2.83148088100129,2.2677790647491278,2.5919030632753217,1.2078440138457052,1.0178382601085763,1.9623730177075309,0.22028578787483821,-0.10425955727055841,-0.20321442615322363,-1.1202453577904345,-1.1320841814635163,-1.0908528423439128,-2.0977059020961284,-1.6934673479486997,-3.6692443390234368,-3.4305968631209973,-4.01899140466954,-3.635969715317628,-4.629778380895754,-4.9951044198140311,-4.3645335090551551,-5.0636168166894144,-4.5304357508454975,-4.2196259233261166,-4.5886235959778849,-3.2264437809374757,-3.0637464296809118,-3.1979458764966684,-3.2152131966766668,-2.2659128812462193,-1.803985900970217,-2.3890684983809507,-0.23027098329094842,-0.25432459561665172,-0.951153903345711,1.2913117270147183,0.57495540842034065,0.23201579455678789,0.78588847730548639,1.620423358839099,0.96713736020847463,1.6050101854303682,1.902663808105066,1.8106519081196062,2.1550102542866716,1.690716595967988,1.882546340261084,2.0365258765145855,0.74206254849376374,0.78725015019252742,1.0963234969738562,1.5485842576387927,1.4900525316267565,0.90993728648145467,1.8306081591805894,1.6564587872896821,0.93754662920288445,2.6486684161054015,0.788867892009893,1.4758490025294326,0.99244975021485649,0.22766657317938832,2.1206067821769743,0.76115824413524757,1.4599876528316358,2.5271692136385941,2.0400576797519236,2.05872036676219,1.9515745960360789,2.2882035096778273,3.2359465909577203,2.1634108977945949,2.3580591458129478,2.2825117492769045,2.3075502069207889,2.3884644088422173,2.7000957700712922,2.3035243407241865,1.6997092649205767,2.6788207772554977,1.1443549834613387,1.0963307933070385,0.94383625896305234,-0.33044992417386865,-0.41661009087956069,-1.0045482443749543,-2.1510417855187165,-0.22111832089214878,-2.1206823469047915,-1.7939155498070747,-2.6456993556384476,-3.5651514339669275,-3.4673787996376588,-3.5816547758605086,-4.6431653127280246,-3.7053549342869552,-3.8786065044996447,-4.634499365983733,-4.32102323643882,-5.2122215409766905,-4.2241245608673825,-4.4151980375973991,-3.3243536962826084,-4.32532146358689,-2.5050035230415242,-3.1167616314920865,-2.4793785936020236,-1.5199700469406188,-1.2117342495352421,-1.7045066960329156,-1.2216917396671265,-0.72237367867137836,0.75904669593659757,0.61927673699341779,0.90090345258376481,1.5639250367327007,2.3967203825979744,2.106836704833948,2.5952914114468157,1.481105204476312,2.44509691523756,2.3270300118655038,2.564246247926139,2.4703747744322939,2.5550210995981524,2.3865451536730626,2.7838482440019678,2.9584588463667365,1.227737133571573,2.6031287553153675,2.0838438398501649,1.1965859708781288,1.6467144483870118,0.69278978793578727,-0.2122024656322945,0.60954512562543328,0.93538246676029335,-0.76383619710947148,-1.0388556040757759,-0.093142518690743636,-0.75952905218782507,-1.4752074969449418,-0.18189584951783333,-0.84649152043772335,-1.101725168232784,-0.898406550813796,-0.73316475690655414,0.12937970168314089,-0.75610732168169725,0.22754919644143096,0.67019827402269261,-0.089254332174813322,0.672032982869575,-0.17203990497017341,1.1333209135344298,0.32955954120508713,1.3712339753316731,1.2306954972523132,0.82035952974349025,0.26338095834066966,1.6113424184568146,0.93003746505568852,1.5683363578661256,0.58078329483746571,-0.49119933094414481,1.0438804005122813,-0.16125215597220902,0.22291634590684359,-0.82711083746118708,-0.38318577942365761,-0.48449026861250355,-0.59157243160923834,-1.0217452667797817,-1.4360900681698827,-1.4539824893700166,-1.4473384516056493,-1.9983125503973904,-1.3402538015998684,-2.0885165889911117,-2.4039703749804371,-1.2181145509250846,-1.0098530386691742,-0.76610468493158135,0.2774607677484987,0.30496672514294287,0.94171701993838619,0.45819181522870761,1.1499375483665961,0.749136667447984],[1.8177587415455023,1.2278320139556813,2.2734865182181681,2.8044346138409653,1.8775632362869172,2.3229858764720115,2.0609220815328881,1.9687018957626905,2.736062776948371,2.1940193009419673,2.3715009593849117,2.0100907202197908,2.0617333213803524,0.22345565251786792,0.57965546253922362,0.72853182913216352,0.70214559887880656,0.81103249341183059,-0.37451874770278282,-0.64403570001777843,-0.76173160579857879,-0.535596038845517,-1.0766284505437298,-0.79282878170698268,-1.9438225893186123,-1.5285500615572265,-1.68923068885042,-2.8249918304826043,-1.9460636507289406,-1.2870717009542039,-2.7907476699535327,-2.5190532316338157,-2.6932180552433644,-1.5826838259013254,-2.3686468523085642,-1.6416901233204446,-2.7537321652600779,-2.7469011606702685,-1.2347908432049512,-1.3314150634405153,-1.3823870911823566,-0.93573862704680932,-1.2872153864827194,-1.15968550387462,-1.5841489700565965,-1.0685147881908832,-1.5600446712655263,-1.4019451579876041,-1.1515386747044336,-0.80854043356391481,-0.048721972510706646,-0.63903635795449676,-0.60566420427741274,-0.51887677882741956,-0.16868617046366985,-0.45581030615459373,-0.68243438547335711,-0.88912249583871528,-0.31989342464199622,-1.3914787095545158,-0.44345865782934479,0.11082468142085528,-0.35206076343442638,0.1927542798193187,0.15267856853075695,0.54017377231142438,0.19999401689235807,1.4975050254168991,1.0890437655516034,1.2009355385885363,0.78228251935286053,1.5318338278414012,1.1559848703565052,1.98005339604449,1.7459854703448214,2.5508436111838173,2.8627241511087584,2.7252792759722042,3.9379414420032837,3.9291168441273712,4.0407102776344814,3.6813417091742124,3.7891435773288382,3.043123235439535,3.3688992701321725,2.6316277733134461,3.1619167294376447,3.4230171541318573,2.0744563069470052,2.25932569717782,1.3391061160651585,2.0573330431981045,2.0299901151247672,1.6301091587512564,0.17818729791921517,-0.03641965057232005,-0.10575596169801624,-0.6191775993820402,-1.0409439399559683,-2.3782896276444694,-1.1148881807541084,-1.8414853624171226,-1.6385164445611797,-2.0803757966939522,-3.89319472570592,-2.7077424915310209,-2.2778795236818681,-2.2224053789099973,-2.9382120961797082,-2.5155663659671377,-2.6164252892382316,-2.1998547915345341,-1.5731509096422236,-2.8771240820071369,-1.2141106151627667,-1.9799866056126896,-2.2948107306548224,-2.0980541030306568,-1.425694773925833,-1.2058552561089244,-0.86258251535289565,-0.45458761847300572,-0.49784827188108205,-1.0151221877430006,-0.68849101764351217,-0.5818133935399078,-0.37646285900100074,0.19667883914445755,0.26873565858472087,0.043151308728340304,-0.30378584073311377,0.053018274369934254,0.51827893204728359,0.14457174259103409,0.90760527046905237,-0.25238330048307772,0.26762768484599203,0.6573848606802376,-0.38678021973871596,-0.01961563843846633,-0.6277136023370844,-0.50608304082344713,-0.77371430012905285,-0.51468560015795184,0.15274158779587443,-0.1207643512279639,0.3355083707598514,0.6164385881218023,0.35541579924275685,0.6450202395358472,0.72760815301766624,1.6290181193778763,0.93396559452878658,-0.08378695951180315,0.9156758101679906,1.2025825949954223,1.4335578324023266,1.4532524972050689,0.8672869110689172,1.1790042630702673,1.3706343457274193,1.3177549261499506,2.2136979960692127,0.36274549317989435,1.5347379808188102,1.6182965906460987,0.6250767144134608,1.0420220944636316,0.60404664432846766,-0.030715646508740424,0.35994657837268768,0.95695572411091478,-0.26987067640330376,0.047959566576531221,-1.0199573540073978,-1.107170126021696,-1.4420820363800313,-1.7010089162857822,-1.5968745359964267,-1.6944924656689579,-1.5585687367387782,-2.223176982891923,-1.6874852122031276,-1.8213616081233486,-1.7969037242273682,-1.4906987837313985,-0.64529558401338849,-1.1004118249080228,-0.72425942091849027,0.048617972166767887,-0.18708131899235481,-0.22969540826399548,0.49951513552404159,1.333502190014189,0.22827436066222456,1.1800450396171793,0.61357715540373126,0.65140992324817981,1.1772476459944903,1.6997845773996219],[-0.66214954700238526,-0.47533626898045639,0.30956163511944196,0.66073965480027841,0.782841740719751,0.6445481241233344,0.37642830816835737,0.7748751663797806,1.8011551137627513,1.587460852365902,0.61122074063108578,1.3855818047096609,0.39922416738503291,0.90795035993828954,0.81799925154351416,1.1293811707246997,0.12963653094721828,1.0250146535690932,0.758652273510654,0.78541984476143822,0.038586879229707272,0.27334485046309842,-0.39873478906525006,-0.10352888768049312,-0.34339923437540942,-0.16552475797596619,-0.0469925222351293,0.14482161353853912,-0.9656937620089644,-0.43849350548613641,-0.195928451059578,-0.40965195538734894,-1.8335030695429451,-0.80869865907511418,-1.3877178689011536,-0.72936191032462017,-1.0583654644739977,-1.7849725436361037,-1.546599758113997,-1.4788307345814222,-2.9992618031436313,-2.0911834155206797,-1.9785085963857736,-2.3311609782520368,-2.8350292181844328,-1.2419243342844126,-1.8962374514292506,-0.916494928642818,-1.8825008107132657,-0.81712626238456332,-0.823837331868989,-0.77119327352401812,-1.441609923743564,-0.88129612642679456,-0.896796166496056,-1.0136146372057855,-0.63202585343493489,-1.1015897946959632,-0.25826281605602563,-0.17510321115375543,-0.39801172212105895,0.42062648810215486,-0.34813138482778705,0.11617630029090845,0.55696861331625036,0.4460310188678675,-0.82293772443518742,0.018273837253115432,0.22372075462831911,0.16725569021603415,-0.21917369958503485,0.44684711221051793,0.39484916455220154,1.1293092484979486,0.55840461360819793,0.870976137976594,0.4731990865933986,-0.34908093448681488,-0.14199565766109068,0.25607643485911175,0.45833134485287719,1.0129458765047588,0.36287439011777306,0.52613221868615545,1.0330644835191067,0.69175858753575259,0.42790483434569737,0.83459093684312946,0.83947985846465245,1.5899674592778579,1.4150865458415347,1.6299823851776598,1.1714574273021108,2.4937567353608978,1.5122544211124616,1.0069865089518002,0.80006269844429911,1.491057658596165,1.2025911875867108,0.67816894092064106,1.0514658361603235,0.38903512638238014,1.4460011767092467,1.9644749459168591,1.2379335121252266,0.57035018272794624,1.2666265291233318,0.57994331979545155,-0.83249741063169147,0.0445149042694606,-0.0043525625868637632,0.078922056744760827,-0.438384636534603,0.083267828075418771,-1.3062268299071875,0.46001608562926122,-0.20839181896787673,-1.1086099809744971,-0.92654196387976562,-1.5758990206858627,-1.6570558189171831,-1.1885835751660105,-0.75178864738382212,-0.7830857955230679,-0.16584583220913718,-1.36538899538263,-1.2715073228990093,-2.4729391952168189,-1.195924860378051,-0.89650853582877033,-0.87274763566577007,-1.1337808717947497,-1.3704625074385985,-0.84965288518906068,0.53611861868142452,-0.95321171384701309,-0.41813780378115795,-0.93681417270910661,0.27672621412036191,0.1301677710759358,-0.15103713986862027,-0.23414175831261339,0.53885118396771436,0.087308176851172759,-0.21453557094216949,-0.16496934550429318,0.50744094212761559,-0.050245525587994166,-0.20416809221285986,-0.098867334010546593,-0.47882433970837479,-0.96947230737450918,-0.3364723571705619,-1.372090139126416,-0.94317689093663926,-0.095174260553951728,-0.990910138814582,-0.51823617496958652,-0.75979905211344911,-0.500779556893452,-0.40308093189147565,-0.82859993530072007,-1.5917211133288336,0.041342689109188591,-0.67723066655471809,-0.17515933724002186,0.2564201697261636,-0.67334667622101774,0.597522641481815,0.44518866083426623,0.24266608852652705,0.10872803566322289,0.99521843229025708,0.71598320727255937,0.85391748182713967,2.0213905229503837,1.0305348551885094,1.4095542324350814,0.85282693123992459,1.8242216842642738,0.90804015696393292,0.47799982757633935,0.97461488265170748,0.12029932474995719,1.5911476889488285,0.34389264394758112,1.0411942219496941,1.3141162866455671,1.0486066156636586,0.939696868723299,0.26352482389967485,0.61001969129097189,1.0452409695925764,0.14317057853916976,-0.143937829374294,0.079012004704302974,0.81675024241133953,0.22122394239321641,0.44750282955758414,-0.54968926292856024],[-0.169370470720412,-0.0087258685205714448,-0.15017781944271091,1.2241768020225572,1.3186784120309127,0.17801893250172707,1.2449184261454214,0.78516892851497033,1.31698099890155,1.5332468822886236,2.4495820353669928,1.8340954662347593,1.3409838812027224,2.0689310730085237,2.037328814060464,1.5759771549972379,2.4971970961309049,1.1250919157641341,0.77684597531435406,1.0461684692705004,-0.13838791737377587,0.024640070915931317,-0.33843161547920253,0.30576294415150551,-0.33687943621664013,-0.69634715217180365,-0.84023173511212512,-1.1339180641861542,-1.8277605034369835,-2.7835746986365137,-3.1132019657864727,-2.41190184178355,-1.7141138940586242,-2.199509143206706,-2.112556314262811,-1.7122347742730621,-2.776127712685136,-2.7278515984432437,-1.6930686069541356,-2.39451746541195,-0.69143961768381823,-1.0703450038191051,-1.5323711719012814,-1.0818433451217435,-0.30310509071825242,0.1275812239160784,0.68716799794709715,0.594039842588119,0.21987772228632385,1.0523446559370961,0.48535385779762874,1.7529302722915241,1.9212636229823725,1.8343456478488935,1.3766310030138293,2.50571477659177,2.3435060509459009,1.4254932541331253,0.99948169192577474,2.010224422178116,1.6140410204948987,1.0944924789413824,0.94167238096584227,0.8770900027212345,0.011311267577266387,0.69331798106127907,-0.22241496478979583,-1.0852698746088567,-1.2420149942953935,-0.51638157422290731,-1.8406597312817041,-2.2610420269208795,-1.331530000279586,-1.8863186452586622,-1.0371075428182073,-2.3195806673003996,-1.866235780652457,-2.147618329136705,-1.8287591638692602,-2.7032959400545966,-1.6070197570985945,-1.1377920161490056,-1.2787943592666422,-1.1871027533487808,-0.70657854376700924,-1.4639436491592488,0.5075452268663494,0.3624981339670218,0.76373120809411521,1.0985746658271611,1.1772311002668974,1.980031149654379,0.91816363701093218,1.6129573633188341,2.3279780559614007,1.8752281401711002,2.2605502878012507,1.1457941232472262,2.1078689311384222,2.4161362193536178,0.66142623319048566,1.5406242007094604,0.8832261621555445,0.82558316451567193,1.453029182261155,1.5085935627833522,0.66253420632380955,0.3143995172801744,1.2906615924293363,0.17629742603722914,0.25881676639632439,-0.23253982700315798,0.095462542992402033,-0.83857043839935341,-0.25367547656961825,-0.71095677850688566,-0.33421110785644431,0.1684349157732698,0.89584749409445,-0.18978487939273675,0.44673433411675911,0.11476918515524975,0.56537918831109768,2.0978078974836443,0.9854923264044313,0.7808571454201374,1.7548102744352287,1.8741238580072332,0.78660542354712626,1.3353232090155647,1.5065362754510172,1.6419233919793994,1.8190691142725397,1.441477745638738,1.8491599957157534,1.7339453726394882,0.89945621595397418,1.2887086124245197,0.28051120448551081,0.95911623403310187,0.11431317079637607,0.23077078035344883,-0.99516388344508444,-0.48848968783655833,-1.4030302417952571,-1.6835999266548083,-2.1547780884970686,-2.0366287312413789,-1.7986297518627516,-2.5472376552898264,-3.2371628746441266,-3.8380702800418582,-3.5858081210718931,-4.2376891475986138,-3.8385138815317581,-3.3232025227954325,-3.7805438331449248,-3.026991532922426,-3.7051097327538725,-3.0191004499492005,-2.8080795921812873,-3.3669013593354009,-2.5859460287858616,-1.4518025521630522,-2.130666417309949,-1.2317059489436235,-1.2290660103868116,-0.81987616736437774,-0.38798782685299393,0.87091152884107292,0.82614340600073877,1.3441660029654077,1.1790324568676089,1.6545850708406924,2.3528351755356582,1.462691712035548,1.8524481320659951,2.8201191918342143,3.1887744901613582,1.8971115419770337,2.4937132632390506,1.5579432995213398,2.3840097829995357,1.9671720812635063,1.21658016818099,1.78437481429428,1.9584648027286091,1.0400585742971105,0.46482890051552817,0.15912799868148331,-0.084553396465225528,0.57327055306675145,0.40515207381064888,-0.26757565662067662,0.42003587859239594,-0.53914700710659325,0.720545224365586,-1.0714014694442362,-0.097635364935665325,0.49174037286182365],[2.3721355777498752,2.6086986113109352,2.5235144383573602,3.9234405144007356,2.8375899456459392,3.5027475419820204,3.8003419708994772,3.0176154776720541,3.8542691266954061,3.2429640228221674,2.1168592149517389,2.4148009786226852,2.3198276386895622,2.4092222350239356,2.1048385591811893,1.7211812518995808,0.24565431313004016,-0.21775101150986587,0.78101910853349521,0.078398540987919263,-0.56521206321671058,-0.42420320555145796,-0.67959969284431987,-2.3668798546178325,-2.3019309108614738,-2.4072831043984677,-2.8144786244456821,-2.9721098525765042,-3.1184745115331709,-2.7548497605993179,-4.6942183786685856,-3.2835266796763252,-3.3066677766744572,-3.2592172347382871,-2.7424164007053915,-2.2877332482200847,-2.3943235652861072,-1.9030198725475933,-2.3678563929755949,-1.7093615759469003,-0.93502085501788146,-0.42506488492089289,-0.46877448583859771,-0.1332893925201778,1.1355350660750003,1.2407934050603122,2.1716265121635718,2.0254155393847633,1.3937294531612285,1.0420204632929102,2.3437220448749732,2.2298429041394296,2.6362114658809084,2.7785474052518766,2.1332860326476717,2.171956581550889,1.5576371989156419,2.3013842450803823,1.4766411153284915,0.90839164123184146,0.65091554357233472,1.4819243838904668,0.51542593033535744,0.16308488883647382,-0.17749842069308056,-0.40129312793781319,-0.652489812155714,-1.117956108476799,-0.46577898853605693,-1.6060820788945482,-2.0347017834841887,-1.9062988896431385,-1.155273545089285,-1.3633460152317063,-1.1376362810416589,-1.7649710584871738,-1.4849505223339603,-2.3425937472062355,-0.81504789806094058,-1.1767767463483243,-0.91871942865507072,-0.25813979380778007,-0.72010808164668827,-0.1927197128354004,-0.60345770540093979,-0.36355410063520105,0.29895409174676307,0.02313748180372148,-0.58064814637072193,-0.02333483695065719,1.7130160069084246,0.48655102274546297,0.29974156216921855,-0.0881390256925732,-0.894191321122472,-0.038788874756183381,-1.7850013239741247,-0.85783446554403886,-1.2804011814778065,-0.8162823329901604,-0.32355667868362326,-0.74086570327147361,-1.3408508080354014,-2.028436645482858,-1.5486204571137043,-1.9077869075071683,-1.0954491535375308,-1.5574461995730358,-0.46796826047115681,-1.1030542725728432,-1.6283115523852238,-0.61156685206322658,-0.035012633947508531,-0.34789461553664824,0.53446142572839994,0.747031047437615,1.1413532827357336,1.1705681185529695,1.6601513866226687,2.0444660975555329,2.4243600393193763,2.1189689281452466,3.1692115350123657,2.8802997424767436,3.7952585339770062,3.3560553317770756,4.7753679019967343,4.6013747442643957,3.9569709062333511,5.2311082432113167,4.6762979994277574,4.377033631633287,2.69517837600341,2.8202678441679225,2.7680472586080458,2.5710426342294834,1.6507020162571842,2.16860932071318,1.085162125431002,0.74856417794429742,0.029806642229219193,-0.27010792352502472,-0.063810123440058431,-2.1740838270185927,-2.2250891263221426,-2.4339693384023935,-2.3367448495909353,-3.4963723363327532,-3.5677018140800163,-3.2195656607055296,-4.3937037943511053,-3.404932216833144,-3.84659909525868,-3.7622491821396968,-4.0265433258547167,-4.4350397733884,-4.6026033667641171,-3.5469064916678992,-3.0465258783113813,-2.7892278549139,-2.9145440998149681,-1.7479736596027795,-2.1676714414905183,-2.1446928912035284,-1.5755032343080768,-1.6295550818166467,-0.68590103887106468,0.020552779064864191,0.0093804488496808114,-0.16777666881989511,-0.78471417369826535,-0.2679503431548127,0.37300826955558775,0.22102003609391382,1.0064322947417292,0.671861448663692,-0.37575200603936532,0.63455213920769782,0.61113540721424831,1.097716771237798,0.20263867716841735,0.22664044044889281,-0.75438524464516787,-0.25030880375357911,-0.80700740377047531,-0.15639883321642428,-0.47248494329221485,-0.70204452478259283,-0.0080403706737231968,-0.34498164969881118,0.819915166050942,1.4394069284048285,0.569324631621323,0.82657713018124712,-0.079426402030916221,1.8804180079691968,1.0645468394982682,1.8238766338859431,1.7809549684333614,2.3014168380103239],[2.4372466277177525,2.6536771463686524,2.1323587485350743,2.7736475975116806,1.9721327785545921,2.636983846025863,2.7238598375053491,2.4019924090082028,2.581652788139476,2.7099156195457894,1.2927152747460149,1.1610820549261027,1.3581486233216644,1.6454427412441428,1.35557531210187,0.74874643722108358,0.95678402280172081,1.4473701185066383,-0.38511806910914592,-0.51297593464989288,-1.236212301685137,-0.4691262311660872,-0.33356849863316584,-1.0893859615828383,-0.50672569752617314,-1.2894208230700637,-1.1463324120682812,-0.82787854841045949,-1.079587198186422,-1.2270545595739413,-0.65670559082230451,-1.787262720571116,-1.50012864857651,-1.0367816396406417,-1.26283361069131,-1.2509253323399798,-0.48187801015280674,-0.96058525983877363,-0.36883207802418139,-0.51576498307905916,-0.50092675602947767,0.082501730058653633,0.12949219418036112,-0.8820523670378394,0.044887360944705051,-0.22817296368066481,0.44965573919370549,0.10257887709077329,0.50646155010844507,1.1771974275292729,1.4231997520024802,0.97933766466872429,-0.22431613087152324,1.2913577506744471,1.068277692853107,0.6222243926179154,1.669213437095467,0.17766510164215377,1.2568133203343974,-0.24187514782302544,0.065510647231508429,-0.36645382146741345,0.46113768645527997,0.32721500116927205,-0.20875228620142122,-0.60834724134656848,-1.2685861772336708,-0.40491603552038413,-0.6476453671908079,0.40308029525510708,0.61822279525566182,-0.91194594922659844,-0.77387813665489924,-0.1296828178325177,0.59087193080622691,0.746752153857268,0.8520009686250688,-0.092948237950840812,-0.24484883851288386,-0.46541650851568128,0.282048768781892,0.33086054179373353,-0.88849033310137426,-0.2790917392008902,-0.57419224460461527,-0.11971398096562696,0.027130168556971146,-0.079301970674655609,-0.57783984886249706,-1.17945057419608,-1.3607775844749548,-1.8614252476631772,-0.8417241059265721,-2.0167824633891271,-2.6756978864789818,-2.2373234646178721,-1.2041663877257247,-1.890110598336274,-2.2647862159929697,-2.0195762280510965,-1.9410136252017791,-1.9011330970138931,-2.165136432938997,-0.85685608870945762,-1.3629223964780488,-2.9926686237725977,-2.0621051531551866,-1.4063283390092054,-0.91864835747937312,-1.0627541398227378,-0.67543991254816427,-0.52951754634618275,-0.72602675502620662,1.3059093074898422,0.45715720914335212,2.0325702703749893,0.736442464703469,1.3560381887520587,1.1206491289269551,1.6208549107917425,2.7524732441803592,2.9077559250970095,2.2930034888199584,2.9768303929068005,2.8193236271806663,3.5574457956562191,2.82794138131913,4.0238904622120213,2.4687707558423209,3.3525117234373591,2.38032054476109,2.6719698594931454,2.3654606653194312,3.5192935504251706,2.4813904761390808,2.5227475196839522,1.4189676069863275,1.4264735972107716,1.837050507641224,1.3622457786650488,1.49198742973008,1.4887796691826547,0.20289925673242598,-0.063640570374648764,0.41786442557885695,-0.091097134233416,-0.23170797788141476,-0.87405068845105849,-1.0834170459942878,-0.54784255176146746,-0.13745962818240365,-1.2215471564288034,-0.61005048095951231,-0.58113526938582194,-0.819021594193803,-0.67436903589162878,-1.3207535823457606,-2.3663951112948838,-2.0234117703303633,-1.7356937007094895,-2.4273081829915797,-2.0036283446556107,-1.1426747407800615,-2.2078476844793888,-1.7533090267628972,-1.7511549592621938,-1.0058808355118614,-2.0904208838444815,-1.6764151792740714,-2.6469360720918318,-2.4983199392020818,-1.9198024363042947,-1.8374136587672139,-2.3435897501981429,-1.9676071220991009,-2.4265013603829111,-1.7740590101179108,-2.3676030088846449,-1.8287111485943432,-1.6974967672963679,-1.0570321865439998,-2.2077511165115729,-1.4052132429258073,-2.2140279004105374,-1.1063836412631558,-0.57789130822348023,-1.2866295785980908,-0.62325724501028035,0.33899509743425471,0.80826090833180664,0.24432811825021725,1.253071334601251,0.69808479479331731,0.60421312854903531,1.2458756369399717,0.67207279684739674,1.3610164390507116,1.6758715778759181,2.8662042433063282,2.4538347482746286],[0.28555959074185677,-0.42282880309970805,-0.36081716887431958,-0.52633356075447224,-0.72333357117345232,-0.044262211639619145,-0.15927023925030609,-0.58552867760605709,0.618399716209174,0.44666889316406022,0.82556939242942018,0.59370936858918189,1.3466953699664699,0.84253553948090676,1.4424026946342343,1.7036778003702611,1.3921664397920768,1.0583325700403254,1.6786090006360284,2.1149263482221383,0.48030579541570551,1.5954874419300726,1.555115513950232,1.5837347345975112,1.1235144558580812,1.3681026912949048,0.79277460838886959,1.0558445242281531,1.3514106542850115,0.57526945727198908,-0.07795540052047778,0.25222166760480058,0.36797683058761582,-0.12899835818397509,0.70329166011032707,-0.46713318960321815,-0.56146380823629727,0.382507581468524,-0.34573555014733554,-0.26826747572556331,0.1125602816769431,-0.87752881220341683,-0.616231974201837,-0.040802726349085217,-0.806235410583252,-0.20296754061166022,-0.80310994589606322,0.14698825608484212,1.1335718882693766,0.06427183319568118,0.077140083547013988,1.240957092707371,0.069443050754714519,0.47064407091290594,0.75224936133386822,1.1680517014575493,0.39946264721274893,0.52074355888238633,1.1968510672263362,0.96802282612436708,0.21963906205302119,0.16762948252189214,0.11207852318554623,0.91439961037208728,1.0362833972990226,1.0290472250387928,0.0330314566520411,-0.587639930046981,0.55574038505526446,0.050558827987097051,-0.08870301765885491,-0.45751592812709863,-0.58363655193040243,0.35422890730886625,-0.83457413227559218,-0.65247451595939432,-2.1302492143867573,-1.909150883018528,-1.6593908298702842,-1.2788483854097261,-1.2347039697003594,-0.20293390297001634,-2.3577094667280711,-2.3436593998417341,-1.6009239509287412,-1.8573149371998898,-1.8874058941884202,-0.82248986174867889,-1.4307765499463438,-0.89194663713724909,-0.92627313830271374,-1.7259615179148615,-0.43903020095882156,-0.84789823201667591,-0.9526922800702371,-0.60657806653651214,-0.25032359789386494,0.17802071661186669,-0.26275949713700508,-1.128390243691157,0.18538893401399753,-0.48433629801322226,0.71891781122478327,0.096262366128013144,0.10197304292506881,-0.42163713000309133,-0.70650387179546981,-0.6890749377600599,-1.0146953402205985,-0.82086883640589547,0.14164470966142362,-0.78363609713141735,-0.541690791391136,0.24973011415392943,-0.78269212946279443,-0.6236879433659358,-0.15045143243415787,-0.35176946830500117,-0.15683683323276051,-0.49487087708373256,-0.29266055617582915,0.72566154698317,1.06374120639577,0.18888284411858133,-0.014158416817475095,0.50555184936960118,0.44879295760880089,0.57006537684174863,0.13588858312214869,0.54010155791117587,1.4872014801558584,1.3560847485270253,2.5059160444614004,1.0682757894344417,1.6593920019299766,2.4828341061160479,2.5490168526656065,2.1507355887836646,1.8971325214453769,1.2118540672528317,2.1283156305250359,2.5399429424011712,2.1544718955179154,1.4121863284488037,2.4302827655188226,2.152548434630384,1.4992875546061788,1.4278203313361837,1.2010635527965912,1.0199069669772407,1.0759335474915535,0.26591779711601249,0.9973704321200495,-0.6080892496075706,-0.4269567152504683,-0.88486973303311511,-1.0207757045571659,-1.6995728141454496,-0.45630524204158418,-1.144586910468155,-1.366353222797771,-1.3507222538910946,-1.8256111209919186,-1.7023132751658279,-1.0922025218177898,-1.5754145239726325,-1.2672359471589356,-1.157139110150355,-0.91669585218903116,-1.6307860731870822,-1.6579114610909564,-1.0100046548707917,-1.5811027046902717,-0.43037460254887183,-1.4647452635913711,-0.29630164908230572,0.050501044684102214,-1.3363051765372025,-0.040823767574970526,-0.32406470657706626,-0.49018272608400859,-0.64469375410688912,-0.60951150844526336,-0.078946889961587385,-0.85940592103549385,-0.86720131651838539,0.94768036270191569,0.15616849664470434,-0.14272726923378065,-0.25634332686334349,-0.7966619492670326,-1.699376778193713,-0.96512682598296917,-0.52077110044515007,-0.58710613145818846,-0.12946536910446282,-0.31498941222512422,-0.67767586110684852,-1.1420401095663513,-0.95154209301609194],[-2.4829356101251778,-2.7714389069221705,-2.5213248939683179,-1.3061707126243622,-1.8126785338316351,-2.233859045651613,-2.2292899948664653,-0.69943907637693625,-0.7572037357827448,-0.42527859161868964,-0.0247773857000263,0.31014172114099192,-0.0014677966291288813,0.67810043098653239,-0.26119239800073379,0.47683646391186496,0.81426875411621469,1.1783366986996309,1.1689020511596144,1.6181656769226174,1.7485214397671336,1.0243969081461479,0.260123026279769,0.79403402651150623,0.99626664692023548,0.77140383412837865,0.097158383381062,-0.83556645943351182,-0.89195976104561536,-0.55086796443787556,-0.20286198754651652,-0.82383634603590483,-1.3463894014280871,0.101595444582395,-0.74439451144656044,0.06130066358419417,-0.10302441735647513,-0.82026722538596852,0.047979273739410078,0.57596172890971231,0.14089094685649611,0.40213715687477575,1.1898072339075352,0.79931897694273579,1.2672845523847023,2.5550211973610675,1.9910464828557339,2.6408781993846997,2.3926685838762207,2.0553790524994864,2.7390470586122708,2.8045828839018849,3.2450136680025028,2.9011558757577331,3.6133127639313272,1.9821958352361519,3.0536279416859529,2.39080673024902,2.8886182300359713,2.1454353999524352,2.5928268374615437,1.2047392216732247,0.85555623531089253,0.83781380773491332,0.70118073097447908,-0.858791466611218,-2.0044028176987307,-1.1862930743725564,-1.3576092234211703,-1.0832849612041839,-2.3030570003595257,-2.4135961177783467,-2.3648821379083858,-3.1996725187106758,-3.5344425848739363,-3.2242833972398661,-4.0237845409588342,-2.8537082361377331,-3.5423785875017715,-3.6694102947949419,-2.683482759318951,-2.2534822774947783,-2.3545492303617928,-1.6877094089382734,-1.2703719574108323,-0.98677808799744637,-1.1055886964268902,-0.544265438722865,-0.63732956350818859,0.69325643925358094,-0.14478457902750147,-0.1130889949735947,0.36411590835966268,0.50587663412868333,0.3481203530625171,0.43978931207936733,1.9893105036834409,1.1602897226044997,0.6239470987576079,1.5379358161546794,0.40659551008543016,1.5016273538333711,0.40981453633128284,-0.029298522660607895,0.7992224814435267,-0.36437467196198264,-0.30738894908199249,-0.16947065324310262,-0.59718592816174632,-1.2529967992573754,-2.1536317054286025,-1.629277603340572,-1.0010462796258104,-0.87099464731100207,-1.31363516944824,-0.95904516648600591,-1.2058985990404261,-0.54519232280204621,-0.38691248251118471,-1.2573155449664151,-0.75261229054354473,-0.23164501193402415,0.36178844987478875,0.062307457268694771,0.0088619703701611763,0.53177141135451633,0.15642192082149808,0.82189062027922244,0.90559184590923236,1.6499611960745375,1.6306356400326569,2.4849410615810008,2.1311564757494503,1.5452819060779768,2.3669482373354622,1.6401087443094065,2.1869837763666773,3.207984408344986,1.6272811461117511,1.769741725313644,2.0309845552776253,0.71122769272493891,1.5840535314552486,0.47850688238908395,0.80812593623391737,-0.16425452661839182,0.6841348995865979,0.14394918760934206,-1.0478178560366824,-0.016490640200980233,-1.2655622732547644,-0.88061988957924875,-1.3122173337443621,-1.4767871462575624,-0.48590366553299225,-2.2506457430712681,-2.1730804504037091,-0.89589068630709268,-2.2280514600311889,-0.54044504594959131,-0.067938492386515281,-1.2054985152113025,-0.19636434775481285,-0.1811395195434839,0.89171353253732644,0.41701829018462228,1.4372289458454679,0.94061005663216013,2.1272466086043642,2.0121554390441778,1.4350971117589086,1.6143213031470707,2.334906872174205,1.385524092810732,2.8491188707614481,3.09355059967233,2.5973076225285832,1.5525169582946692,2.05668851258045,2.8194662516247773,2.0520638632984292,1.3826037761989463,1.2033746929510829,1.07171793154567,0.11429516702057162,-0.42437532174746367,-0.89500876636078819,-1.0600113262311812,-1.6214545585509397,-2.2284049706280129,-2.6218659838515168,-2.1273263075890365,-2.9462961517962851,-3.3207602609623357,-3.6729720298423807,-3.7598050558912024,-2.7630785188860361,-3.2696505209188822,-2.7086253020489988,-2.2927941298287231],[-1.3211079534616061,-0.74163059108125351,0.42645325637074016,-1.1117846094919026,-0.652893667000725,-0.976217795929885,-0.22957582101692653,-0.97344326155058636,-0.99006567669080592,-0.48236070623698751,-0.88173773754740337,-1.2943812997432904,-2.1394992006943423,-1.8846371617829,-1.5931547667604993,-1.0844768686331467,-1.2941705879260028,-1.5396924977425464,-1.4667159815259754,-2.2411032135984885,-2.0984999514816907,-2.4493161646486317,-2.2984570777984814,-2.6907720258242853,-2.1220294033176721,-3.420908554780774,-1.4185713128196005,-2.6493811830288156,-2.998663297297409,-2.7765254965763515,-1.9995575781239141,-1.0579965179300066,-1.8682029035140002,-0.83602708740288023,-1.0733788526939343,0.23040394021078048,-0.086996448811601318,0.24110310896313575,2.1868632780371593,1.697147400859907,2.6621966313067777,2.921428808621521,2.9584469095600348,3.9091495834662249,3.5160235283117331,4.1418177242316885,3.3371707728106279,3.693253099617452,4.8183249848375871,4.3289526189458343,4.3087088053962646,4.6510955537126089,3.0281171488314165,3.9869304902620795,3.39577562556388,3.8891760693222115,3.3396955447137269,2.1902226275041623,2.505726822762544,2.0165920129986654,1.9485540635850223,0.92137601242343359,0.10641304698803766,0.297401204155739,0.15465616833808421,-0.324729168787846,-0.55548943697865472,-1.7003840734598643,-1.964874701227671,-2.1371112032499529,-1.7113640677312107,-2.0476116766140566,-2.1346666640070158,-1.8737958381709667,-2.0692392353188733,-2.4762578203076608,-2.2484630992818251,-1.1185058719923449,-1.3830940274854482,-2.1731671785693254,-2.7056530260527056,-1.7796111051038113,-1.3182135022965566,-2.1841209410713907,-0.33784213640952743,-0.093187858416167391,-0.63691428626209878,-0.67451917576507658,0.027111496668881685,0.12109517363450462,-0.23202893074419534,-0.39845010140791687,0.027317466890437597,0.81998261997744559,-0.87397148466061625,0.65297515338375489,-1.2043556948584093,-0.80159648603118,-0.20501289671434797,-0.0539845677641706,-1.1592171272514089,-0.93511819281872877,-1.3588644776342282,-1.8157416678795064,-0.55114293123849123,-1.4167237534941239,-1.1894923357682008,-0.14843290678210352,-1.003538363378305,-1.1250838003314609,-0.65715004231207708,-0.77079405483649033,-1.3377238931868789,-0.18329513777650247,-0.33780278377127165,-0.1578166598421446,0.51179794838952852,0.48287928330391272,0.36974798824225941,0.2968549146007885,1.4749545916311901,1.7776238431488312,1.6997203685009012,1.8320873289668029,1.8981162592366108,2.0045995298239476,2.0395871370669538,2.8993154916708677,2.7556780843621418,1.9356599165744386,2.3924983114463387,1.5881851540056495,1.332119858877377,1.4020210763603584,1.3037306947797054,0.623659180434988,-0.13173167224750626,0.65122001648977879,-0.3201394375458988,-0.43762436151298983,-0.50507428046734437,-1.0927657158321333,-1.4413501150114896,-1.5661899206864949,-1.7301839031031245,-2.513694203986129,-2.2263780817141248,-2.2276162003282116,-2.4315024904053968,-3.3921504412917782,-2.0249206734945067,-3.0315518101281871,-1.39771800184083,-2.5085557051239218,-2.0471697237177091,-1.108819723845579,-2.0106611575872289,-0.74605949785744929,-0.23164811609967406,0.06531886446205043,-0.1617949911975704,0.8363708300288567,0.51967852741723009,1.5453487945075184,1.4849330586472953,2.1917736300919333,1.990691805219108,2.8662398154918951,2.6298391112575592,2.6174675827758396,2.0566789035158939,3.4591967732276321,4.1043578860225445,2.9166669282523561,3.8229293558940975,2.873141421442432,2.3780703526811462,2.6055993117670373,2.10703747424584,2.3780287929734643,1.2999859419131037,1.4789232651994064,1.2352388027827628,0.8313136560440465,0.59295001291591465,-0.068397242151833793,-0.88510806811802178,-0.24343699352049217,-1.1258057124727499,-0.30705070191556239,-1.1966738001293225,-0.75657369471318814,-0.43702597908279373,-1.5464272798801268,-1.7229551346564829,-0.88144398196833074,-1.2669854476647981,-1.3714828655481563,-1.2644735163147234,-1.6359695647611545],[0.90926572575137055,-0.023063904271610036,0.57386307734395769,-0.42203599766874444,0.34270341164809792,-0.45009703737309165,-2.36885844728634,-1.1311894453562381,-0.63372448611513021,-1.6460880250170742,-2.031596251794046,-1.6955806372264324,-1.3931675863917459,-1.9481294071612139,-2.4424547081418946,-2.4049539623597909,-2.54781519597086,-3.0200750645338452,-2.8306797264798473,-3.1099876117587737,-2.7997862284144008,-2.8586535006967768,-3.8724108059455435,-3.3917678817364956,-3.6633039853130041,-2.9001293818803537,-3.9581870112905566,-2.7216578339505682,-1.950514616018397,-2.0539183209970191,-1.3365188698965973,-1.8000820721312603,-1.795263192739891,-2.1634829828372117,0.36255607723951022,-0.511664133109062,0.17774131305938656,0.21492898378254532,0.052980478591314745,0.37980841581662894,-0.43920122437991982,0.59081812934561173,1.1374963926311397,1.2568255372149266,2.1691402216443763,2.2298352483959518,1.9213861080464125,1.7307606159585074,2.8375293334984737,2.5045807263112065,2.7845813651184561,2.8381941480396633,2.628727392407149,3.0464757786008168,2.7800861494248963,2.6464905780208143,2.0287986804229905,2.5018719466808239,2.8458057233978851,2.0750487210435313,2.3606200489097668,2.3487798031982074,3.5602437221569083,1.3936355510786846,1.5047277958934315,1.2430691701455725,0.54169151771562229,1.250204501746464,1.1722784987021106,0.68141024824305,0.69572356057093165,0.52552971135576532,0.28740094929116389,-0.12307532108608132,-0.18722584681900056,-0.41395370973276707,0.26949120469004928,-0.21424299325495716,-0.15558465703477803,-0.78554549742296909,-1.6307849997720254,-0.63254236145789722,-0.45133223525229449,-0.296540668963604,-1.8797825114060402,-0.917411237268704,-1.0792296617347661,-1.2231125759876078,-1.3834776593234492,-0.431237802459042,-2.0503750247426553,-1.4299165397594553,-1.3592754658741439,-2.6073437439628062,-1.4266027842037305,-0.93122692661482354,-1.8042028724879224,-1.7924992538887219,-2.3459567677086759,-1.4217887093543524,-0.99475451482945132,-2.0657236304616884,-2.0212480527583123,-1.5393266514959494,-1.0705466766145411,-0.42462488224951511,-0.55869593323065181,-0.72549887485892683,-0.94257161750485041,-0.3205715640728849,0.092101830198352053,-0.36300020491467855,0.28714877046901383,0.47793089744289013,1.1707658618976904,1.4414500646918584,-0.12346969793735785,1.4866289650002451,0.73022533653208432,0.49669237097268681,0.92827622269108478,1.050622627157044,1.4815207870864988,1.2309012271926698,0.6533016440809154,1.1229515130145729,0.28872781436400852,1.9571963034704896,0.74741087728473365,1.301248132682487,1.1544396308460709,0.34593276391179284,-0.25062168360541459,-0.42278330210597259,0.56944478148609867,-0.0751111908553719,0.33378095354705606,-0.612214172475713,-2.0177602958288259,-1.2649615557865999,-0.48255454124768815,-0.92176040727457442,-0.22055994088007602,-1.0958941272660145,-2.2602872416964841,-0.19695462196019564,-0.95537274847565457,-2.2537422981345396,-1.4232826359922344,-1.5442412106616359,-1.8676936339117265,-2.3418661909420662,-1.8352087319858357,-1.1330149055667174,-1.4402324206019861,-1.4680730888860722,-0.87937093013023915,0.54520384985186054,-0.48006545924343941,-0.057085339431175053,-0.5471644961394091,-0.88788991468387213,-0.50547439907798253,0.494504339533474,0.409007246665179,0.51758461378983434,1.121119895568504,1.8728176200255033,1.4464084912613595,0.92059412347143843,1.3514191685676602,1.7737434613484404,1.4300234163837235,2.2125094436677593,2.1642801080413396,1.9060990854171334,0.8787214729879198,1.0998667695860653,1.2254785737102551,2.2773418232251323,1.7561385469857855,2.7276237657350961,1.7175699748022042,1.7661606462609576,1.6461789363511654,2.3801386726189921,2.238122308055857,1.7704252696477505,1.3855074389548478,1.3986399910728124,1.3472268561991974,1.8029343871680366,1.8934548596712655,2.1232272204880021,1.7235044468377798,0.55803424601563856,1.4811835651909322,0.91077404921343674,0.35041072044548616,0.45010171512295583],[-0.47724590025859848,-0.50045443478917062,-0.50477597165618193,-0.69588723676158182,-0.88232225743833548,-0.056599413940586651,0.054700020573602681,-0.37797706034920808,-0.58595919963919907,0.12798508940053621,-0.72498542560177459,-0.64749750869644518,-0.42684620247858773,0.078879196061450885,-0.44267673497413246,-0.10111608448892928,-0.42360386396815131,-0.27560357936383589,-1.4068504077085278,0.23180965526320185,-0.03019576806493296,-1.14831991197957,-1.9571949949173146,-0.00808475193577285,-0.93669981952107328,-0.63933624756476415,-0.95698132491111054,-1.1973926324627935,-0.96794509331338252,-1.4630538401922986,-2.6914477371890135,-1.6245400364778391,-2.2434658997376276,-2.0649821070269132,-2.1170723644688576,-2.0934096698114972,-2.5527121567299704,-2.5371419118875917,-2.73346601979313,-1.5595444784677885,-1.9395756972570246,-1.6059830974956022,-1.1302949273036964,-2.0121110743571453,-1.7619789334502298,-1.020803112298327,-0.30699977238316262,-0.79242275295394649,-0.43019808342697535,0.52031782426618212,0.71694844377049494,1.1386070085593509,1.4950777546913117,1.94338748319738,0.80395948380311244,1.9705124483482432,2.2790453837912694,1.4086591226004124,2.4356235425446138,3.0158034157977593,1.7770833953899896,2.5760778840368195,2.2469026901100171,2.8321511556309735,1.8269780511491063,1.7006662457456319,2.7626766578475683,2.1811988866929775,2.4294263682107622,1.9664697187341145,1.6180897817847886,1.141784487824514,0.19102009885531968,0.62234594026527168,0.51101424019551323,0.45599666699849745,-0.44222039413625192,-0.10912661890837533,0.55830031483390652,-0.69749085740893857,-0.4631985755909519,-0.49454832590853132,-1.0816900640489218,-0.45447811819184947,-0.51662960487136111,-1.7739718269188562,0.580915026983911,0.23760031009233973,0.24750528327137733,0.1744276850435425,-0.31435807873006794,-0.62815852669333372,-0.14940330725727816,0.4463322435658269,1.2516260173430831,-0.11532321331909645,1.3939728368771773,-0.009124088972783273,0.955942449862429,0.53354393178288528,-0.41265600585005013,1.7840351387328726,-0.12050185819056503,0.056307482930223696,0.09200846696421594,0.019569835206382585,0.18848637666407958,0.10948754182331771,0.84512094580137254,0.35932704666133,-0.021641544509571342,-0.98078639760041786,-0.25706266810074951,-1.0252841621938364,-1.1879266897384357,-0.87126071886621159,-1.8226905980729682,-2.4863077533776505,-1.525531936036624,-1.7505907982691524,-0.8870896095233034,-1.4047856884614258,-0.48918904031581723,-1.3840580567544172,-0.7004675107691003,-0.64227024207996075,-0.55022143623089437,-1.0782681068447786,-0.4032311154557584,-0.33268786702423425,-0.39467682414384841,-0.10019730602597629,-0.039710162613377342,0.18343680119768518,0.027928032970891303,-0.050982177647558191,0.314031979780359,0.54514246147049128,1.2164101199693267,0.42123003236727463,1.4653551068385815,-0.050020141928116435,1.4985727779652112,0.53092859063233644,-0.76980683177093745,0.519494311127814,-1.1094769789248538,-0.000494887232598773,0.77314346652961174,-0.619442827042858,-0.47456050331800692,-0.55633531364469191,-1.2434314666452129,-0.79276040783001855,-1.5172446465970588,-0.99335529936425526,-2.0071260315418806,-1.5184463178003491,-1.7124972856451448,-1.9365934981515021,-1.4493718472720865,-2.1428842541801858,-2.1054022250437803,-1.8602928728743369,-1.3545106888748886,-0.324715429857693,-0.49185156808577796,0.35476433838968069,0.48570815522428984,0.021971409506099909,-0.1435418842257572,0.823862193304236,0.90295006694194535,1.9444822834486868,1.2051274656481961,2.58082586121398,2.1128821697745712,2.384825923043989,1.9313298039291098,1.8843090722375069,2.373293569344844,2.5059793723312667,1.7638414457693816,2.514129613885796,2.006261692534768,2.5020925822667581,2.6521967621297238,1.0662834850590179,1.4752896016410231,1.51161666901033,0.31736273583697727,-0.066564506858799222,0.16691652181536365,1.2447430342397015,1.3941135373371445,0.061973421166328924,-0.036829625657770287,0.10089392266485263,0.2600197291720231,0.15318670217930358],[0.04432302841852187,-0.75730462704415757,-1.3338842174239813,0.675358232641915,-0.70208662650978071,-0.26754309297756507,0.82621190699605407,0.7399577326751583,0.95028792822814412,1.1604675276173522,0.69996408726616743,0.9044446074652166,0.84665451815613646,2.1570770129814387,1.379208924344375,0.95984446811800839,2.0172226270192044,1.5562461674322465,0.89040381480747166,1.5034249159304443,1.3561410570764605,1.1328066602973925,0.29356645847041352,0.11836789531080391,-0.23764842832724728,-0.2625670131783831,-0.41149225511681403,-0.15851579952561223,-1.5845576852641394,-1.8654985196011067,-1.8843821321046619,-1.4256199753550844,-2.3761250490742074,-1.8401333790683585,-2.7536650233170983,-2.8479878881176588,-1.4998149710700754,-1.9156319776273318,-1.9563514835258946,-2.0894796695021847,-2.2073873067276155,-1.1017156882395489,-0.87845647248397885,-2.059119402818173,-1.4743643926685177,-0.79626651528171755,0.41558658431740003,0.14678233537207824,0.867874014127253,0.16681412517596322,0.90789030210183153,0.53987880425688428,0.057393064118055337,1.2528535781238888,1.1798703651678095,0.79712677723773906,0.87691375703905583,0.69348232997829107,1.2345210928866646,0.30812652625721848,0.94623828001642318,0.29909957810472021,0.10026400329436837,-0.46109584880108156,-0.68961032384319232,-0.856270946118806,-1.220931286786568,-1.3409725661962779,-1.5138581774573954,-2.1296409520354151,-2.3854006870276692,-3.1809988754969809,-1.4043023272602664,-2.178759519662949,-1.546804898670582,-1.7712011705039852,-1.0232322992069618,-1.0734241658452044,-0.30002525434652849,-0.65709802308541654,-0.20769475140874377,-0.33809394201982174,-0.25484448329037829,-0.29295903519823197,0.25783008329978785,2.1977018938020123,1.0334830162804591,2.2872305463758655,2.3895260697003891,2.9925545349058695,2.7900183501468176,2.974892644371117,3.5875241342742559,3.2219356961647074,3.553323462261631,4.1360655041808494,3.675048253424658,2.9497981603959613,2.9368290841831435,3.5605553160783958,1.6914889395045081,2.4339740322280283,2.7054653872563517,1.5682141581465168,1.2732199955767685,1.257064361247801,1.5250723007631772,0.45120878019300925,-0.6937516516798905,-0.771887501044771,-1.4539381357167498,-1.4079996052648489,-2.8962556103609547,-2.3996793827406515,-1.7897352657799961,-1.9548692093715094,-2.9338434820850954,-2.0374158661342183,-2.1702227288044318,-1.6799543370226813,-1.9926584025847762,-2.4171401043264944,-2.2159038411324756,-2.1825659284172341,-1.2222235018858161,-0.68468580873325935,-1.0817500124487258,-0.79511171373227718,-0.68159293068171278,-1.2977171320681262,-0.64161920224291491,-0.058466149503363485,0.44075218052837983,-0.05283858664600527,-0.18428904384868161,-0.11628720758917765,0.80919516881699449,0.48449991323245634,0.31548891893440079,0.49776251221933315,0.37012271509752542,-0.56267354960941951,-1.3219858233183839,-1.4133099540240961,-0.46426556702126709,-0.90901567403942252,-2.063348165241889,-1.7814003454439897,-1.7923040509096142,-1.8599207789974919,-2.0749667339697448,-2.5472968103436644,-1.6253835329432813,-2.4274549983345159,-2.6057414956983092,-3.0175685111747668,-2.5335527765520012,-0.3307406748034325,-1.7782126912416645,-0.46955683054476338,-1.1694496867569462,-0.06765909819851279,-0.52370865576288439,-0.2153907973069345,0.691030173159671,0.96629508693666377,1.5977161469822154,1.3457066920290839,2.2006555647294284,1.8640175969495745,2.88134562203788,3.0129595054630305,2.9944887100562947,2.431345824853127,2.1954041775028683,2.7262956916762375,3.0400415384033765,2.1079370362571077,2.9807150834055243,3.3311332279530488,2.048150682480347,1.9326787211802949,1.18255182121667,1.33669848012345,0.578618485160233,1.1827588587398208,0.69796303080275757,0.085280031945660861,-0.49993880723978579,0.15270628621947457,-2.01481220420584,-0.996455546761728,-1.0320850674544073,-1.4950015618099237,-0.9581166218182936,-2.0625649773220731,-1.6420949357243226,-0.9418685698448509,-1.9764724190571137,-0.75868067286788088],[-1.458532275239202,-0.50227385125427215,0.516506101105142,0.26564102608247769,0.16813412356957891,0.46771192159167274,0.03703324998099583,-0.12706876916397869,-0.31466866069343613,1.1128608960845026,-0.096255975284305517,1.0874351550505095,-0.017490926303539212,1.1097724476499731,0.94096715097931238,0.80808759127538266,0.43611096235540181,1.140184205417702,0.80668252410774555,-0.368410274616099,-0.14951737569421669,0.18142409200329457,-0.34686218357541321,0.25944442758557962,-1.0134197742128443,-0.2090110837717159,-0.17486910884547435,-0.40616376903719692,-0.47079596063787582,-1.261292850044941,0.28785047924710139,-0.5484433951541049,-1.0204550216910206,-0.76790994683997227,-0.86129403226984025,-0.45453801761544771,0.48226116915658274,0.1831559399468983,0.61786538112621092,1.0984244929571132,0.52437383676632765,1.3868603790127998,1.0468822200221792,1.0995591594686744,1.4928563641512678,1.2887970877376882,1.0991349617339119,2.2139831138748063,1.3241744534601327,1.2070426104679477,0.85075765335644282,0.13644302107276784,1.278987319382684,-0.079052041896624459,0.17908878478832096,-0.15369932676579018,0.068864712294311015,-0.96749778629811423,-1.592347163480097,-1.80762732460079,-1.790534600686974,-2.2743018182291066,-3.2133699331987993,-2.2314817891758403,-3.2250144017828806,-1.9222455901972932,-3.3446108580524143,-3.5702214783981696,-3.45506388762586,-3.3892375329792328,-3.174113787244099,-2.0286767955703833,-4.1551173255327924,-3.4808149443745724,-3.060500117098639,-2.2651614643583722,-2.2695272213373254,-1.4306172669561517,-2.0110785493269958,-0.32316276602165772,-0.34073972519097218,0.12398687383341317,-0.22611442154164191,0.51133047068320026,1.2947925106048459,1.2106399209001819,2.3530705115031214,1.9453715772381515,2.1481390556273285,3.0135374895403961,2.5250081620358014,3.6044672748502022,2.4626227159361673,2.906709470885779,4.2985668481612063,2.8220645650945553,2.9388302495218328,3.2001338830960058,3.0965555523473096,4.2130113355016672,3.2570433773498175,2.0957873075953337,2.462792728970729,2.617344810362741,1.6580895325611353,0.24977429735001277,1.1008665836504736,0.40901360758397753,1.1719833232995287,0.66602159079903733,1.1067722417825867,-0.040527280460445514,0.95089495133724611,-0.27731646819807332,1.0108633848350519,0.3137817110279979,0.80644874959954238,-0.50462225466213206,0.64940514960568607,-0.065194829693082323,-0.52168915683764161,0.891542985468463,0.47189655604858294,0.18248477174449235,0.355938208609637,-0.037741965295294777,-0.534554419641615,-0.99040818451848145,-0.89736424883244947,-0.40327682818718524,0.013800027213257349,-0.55188039903201591,-0.63093154867600831,-1.1430085754455761,-0.99266890336846636,-0.98680275649463178,-2.0184088487165672,-2.5334273077374059,-1.8784093298022575,-2.7137987957191565,-3.7255833050952272,-2.7015744209315078,-2.7062717238624963,-2.8032956386029655,-3.0835604788777382,-3.3829976033268472,-3.3131437639084864,-2.3932001220781935,-3.276748198267744,-3.3668844703478387,-2.33897262274364,-2.2539020518603832,-1.7321716346802787,-1.8165513734172485,-1.696459191142182,-0.92555333881374291,-1.322487492441953,-0.9841725854419191,-0.82798017301544458,0.86145595327749191,1.5867302489332242,1.3305038630015951,0.92061887385298535,0.98462067455154823,1.4768760301406083,2.2628909171048925,1.8572939834427609,1.9465330215238477,2.1417163787763993,3.2432997089358224,2.0692752540431756,3.353680843846238,3.2860275681431803,2.4631115610480254,2.3874247500584187,2.3966473536306885,2.1537271028311533,2.6485722478092968,1.8777346809348725,1.2263151589427148,1.5917302465661138,1.5088372674694919,-0.16630080260561397,0.98173907597913035,-0.44013980487464244,-1.184791254846238,-0.43868975413569444,-0.14615901790373242,-0.30289855375322372,-1.2791013057294993,-0.8675630856767782,-1.36502546640344,-1.4499917450256585,-0.79069192149295064,-1.2208365505644545,-0.79829641062156365,-0.94096032209394276,-1.3866146080547215,-0.643102139119376,-0.83092272009966939],[-1.0317271714353953,-0.40405064261615015,-0.58942378766603509,-0.91887439559919992,-0.65965931300052028,-1.1948376411478292,0.13808600796038295,-1.6422976339644648,-1.6351672891520392,-0.34834531967058946,-1.255091282270266,-1.1768340056640654,-1.1500645316002736,-0.34481629051257506,-0.68496140084650292,-0.43554392375382878,-1.5431295005804881,-1.4229996721171316,-0.33201753886314367,-1.2822158071734848,0.71544956213081357,0.03959715798049035,0.65267088451408739,-0.76216829190715163,0.51799873520703366,0.6269717570164236,1.0849519975096054,0.46804488054418741,0.085285349012290457,0.53380838648792284,0.74767975917430185,0.96007560969833639,1.1181088649140218,1.4446725738394748,2.4508036885760669,1.0507309019626767,1.1909893818492714,2.4429706316805579,2.8255649178420561,1.2353000139180437,0.63766601341859452,1.8183913124392594,1.1768245951755478,1.4314195980622368,1.3440861674292865,1.4624103314742032,0.53798120035328489,1.1863059214527065,0.72011339626252591,0.03166176922629127,0.54132904572313234,-0.0855875612180299,1.6855649657095313,0.087397833176624265,-0.18373435514593334,-0.44989783545796092,-0.595843721603881,-0.1303267544121105,-1.6318984747022862,-1.5180984494564589,-0.81719754300366232,-1.2702910728809556,-1.7737897981186503,-1.1610630148564809,-2.0353470549443693,-1.6440396839217186,-2.2376924440977044,-1.4417740378677197,-2.0052458477829647,-2.1581599322710123,-2.9569771199113344,-2.7746469573201407,-2.3675545818918331,-2.400907497921799,-3.2243645440234348,-2.3186951854533464,-2.2071393331701197,-3.1754635893957448,-2.7583731192236032,-1.3256507132488555,-2.1860765607894184,-2.32259380227872,-1.7663313214097696,-0.77904591100711507,-0.99502299142879691,-2.131424426354311,-1.0707079130055641,-1.3731939141350773,-0.66086097519817311,-1.07972432788354,0.43189256450557389,0.838905934414221,0.97339953934965551,1.2734838134838222,1.6142110686303557,0.79710646421076492,1.1550068997208598,0.91281659065747123,0.80899000124309606,1.4993109594337539,1.472888545973893,2.3049782739429729,2.2249875769667335,2.7490759375631639,2.1662891741604327,2.8879507816303813,2.8827378796906666,2.0136846342226393,1.8837193541121713,1.706721583615241,2.5369172779130174,2.8776015387953113,3.0432272749640643,2.2226144359680751,2.4209238103255037,1.3926979650787694,3.2961302392841159,1.8709054348348662,2.2951730002277535,1.1340193125219571,1.5030204954553459,2.0299686755139676,0.73359666624248254,1.368389193242415,0.29966267980506822,0.97858848966657319,0.6617943971226139,0.11046041625437725,0.82568073452903334,-1.1172627882191868,-0.28519595158226552,0.10285216814649978,-0.859350175793389,-0.87897610145829552,-1.5333442328489479,-1.194597513952153,-2.7143282865594269,-2.5469500485650927,-2.6740592607580469,-2.8312569587580012,-1.1721858973321349,-1.5937002369165614,-2.8379967191484115,-2.3552050771946047,-3.4293607657012557,-2.960339481460442,-2.7370160866112463,-2.4007054957213176,-1.8937019803054249,-2.1097667129723541,-2.3324442451810152,-1.5881623886033456,-2.1722626746425715,-1.6090702496168867,-0.76280762584244877,-2.1971165456122388,-1.5752108897753181,-1.7398173697538843,-1.8439309253216847,-0.37033530726449237,-0.900088317664596,-1.3241858163783677,-0.79126011356135617,0.40398814753585016,0.41974613811111117,0.330851322425951,-0.34828552990129957,0.34410026270610211,0.27219713440232074,1.559741559639886,0.68325348060003654,1.5504651661842268,2.0687395594930122,1.9710820653813681,1.6798127662970588,0.77548129728481541,2.6186050170451076,0.817169887040455,1.1416470657414854,1.4692191400288552,1.8297330796893951,1.5591406128405036,0.974672091947081,0.82123980237204175,0.499517725120474,0.98206402325630737,0.94240853043497474,0.922334138366681,0.70169429759255075,1.0382861027175849,0.6064328195832045,0.27028562270628465,0.31004273812383309,1.1522217815119311,0.24899377686764732,0.63237100570124249,-0.20092609255443419,-0.47079096418648331,0.47137623696742004,-0.019240353644348007],[-1.6445815845166154,-0.66318859238794947,-1.4418709974210331,-1.2948434095368526,-1.1920595711244317,-1.1194475154030537,-0.74826366454627469,-1.4306214554024845,-0.24733169882623451,-1.543231588093972,-0.10760692747416575,-0.44182984641553125,0.5660235929426034,0.04289478594021523,1.1805771370850024,0.577689203258744,0.61978567006574747,1.36576944882278,1.1593691524735306,1.1617075431640451,1.2010011997580172,1.9251847913569895,1.2483095177381605,1.4021287098231099,1.0975745387048033,1.219417692233745,0.80416654926379971,1.328872518558712,1.0621439008659719,1.2029196729250873,0.75061225544958365,0.63193182465318709,0.41285083316377225,0.94922145955005344,-0.016799957004402816,-0.24514154077621708,-0.017185865921186455,-0.42656186976810506,-0.93279895693894843,-0.586915571565253,-0.56016289090282778,-0.012830482052505365,0.077456897797193724,-0.21094750893608616,0.11197993725456441,-0.33040156805160176,0.19899229191133125,0.23741975511277874,0.88297720449810324,1.1148186591741336,0.54901301529913482,-0.011467024191258335,1.0986713828243408,1.7474245681026019,1.7906021322388814,1.5893699058486863,1.7832240092571547,0.79921954255661243,0.57899007223078935,1.3436350291189734,0.907659403619386,0.22750518632226202,0.97583317019119731,1.0849841187211338,0.27011153235008789,0.0017357843370524595,-0.37241522199955907,-1.0971588256555889,-1.049068005668945,-0.94839587029282724,-2.2589159483473797,-1.5105599668012164,-1.6415040922247655,-2.5712933983497184,-2.8151078675868888,-2.2677097456661008,-2.8733811106354952,-3.0032728147726511,-3.3709531652982152,-3.7281139695697005,-3.5221229611369016,-3.0165270275879226,-3.566437257833269,-2.9605566643916066,-3.1088196230606124,-3.16853696066911,-1.4063183748886572,-1.4808075164751529,-1.689473872829014,-0.66598343282084926,-0.84002088751042314,-0.28249220490501792,0.057593373184300856,0.347836514357095,0.54355220037000984,0.85251510984005385,0.79313831830191051,2.0121857285657843,1.5098409541429103,1.838656614407534,1.4862675548877013,2.1635718224011069,2.984222245664915,1.7258358354113228,2.6538793181229874,1.537508074360121,2.4827403344050563,2.9004033370119426,2.0859514853407624,1.31923768825264,2.61420973667342,1.0777745696411483,1.1093264355975716,1.5077581079792537,1.3236420986005042,0.95398196522955059,0.98782624436947453,0.96492841312924682,0.56541484573424639,0.70703640580056315,0.33524836719229489,0.2292023239411996,0.2248267957974377,0.53201361640671352,0.01927633735493095,0.22386520969172552,0.57587693397967077,0.88861141136293775,0.68798478931920115,0.8721155812702357,0.63980370659037711,0.26924759340164894,0.79288867467550528,0.84882130955418122,1.8921304671111976,1.6993542574443143,1.4277307942498243,1.743225622000399,1.1035299692582818,0.67970527472851483,1.2659204583899579,0.341952349055495,0.10852044702919689,0.5635105188053684,-0.22627678388356004,-0.70170393302233269,-0.83250501314730441,-1.0244105976942732,-1.6613664191683875,-1.6145763988560331,-1.2400387873549554,-0.66792557575174971,-2.2723627639338257,-1.8841064088286474,-1.4697745563656486,-2.5116739446977525,-2.8014893908078462,-2.576010236943155,-2.60552536080632,-3.1073280664479648,-2.583573821954805,-2.8396326286381317,-2.9547653836793182,-2.1722050989548509,-1.572872228413345,-1.4841425871949574,-1.5170651786934948,-0.88102612478335873,-0.39251157138483561,-0.69341786679503237,-0.47529111682183856,0.35237319092744318,0.29043822726661084,1.5112375371001296,0.18985940588685002,0.32559836458509639,0.74449991311559083,1.5141662633410058,1.5469437910652328,1.5483772956719379,2.1858517789927525,1.6506922137644862,1.9534304137639369,2.4383463836050758,2.096225994872488,1.1323256946175793,2.315512794910028,1.4960284526821952,0.91038649146721207,0.032317145935160885,1.177656854507513,0.02918304392151469,-0.21001740395154211,-0.36756916545422136,-0.48404681744003819,0.31396209214819926,-0.28758937308161958,-1.1475841819869042,-0.42518339205394651,-1.47478264675771],[-2.2518486151329897,-1.6025044216816096,-1.306938637729532,-0.67436654646854088,-0.24273369264542466,0.61244253753593236,0.25044119056470027,1.494341710896312,1.2121690897419704,1.9648224931762719,2.9374859983719035,2.3550678578048752,2.7615107182459249,2.3381719082381203,2.9234067624848743,2.4816519622924473,2.8046613307794037,2.4979411423621389,2.8070892377012551,2.1836482531565058,2.8712991929501452,2.3549519027524366,1.3297638751224867,0.714091230368955,0.6459280055306118,0.62549240856281263,0.77808375902142746,0.67154829674395633,-0.416422398964411,-0.97966821846219754,-0.5782931052133955,-1.3979558772122833,-0.84093180613425766,-1.4367511926291945,-0.73057903391609424,-1.3222431783544841,-0.974307383101712,-1.4505799501984682,-0.018983813630883084,-0.72272992288998472,-0.748285671591312,-1.0872377850746275,-0.0438092326536762,-0.5332779820922926,-0.33122481836328294,0.1412531983124678,0.40088372588182736,0.65005296438945581,2.0010792195649754,1.6800104851987296,2.1887843138361713,0.29025177141913772,1.5706857873310116,1.6701299217975374,1.5743421021992736,2.5044027160760023,1.3195724501704251,1.0264017509527916,1.4014245672509411,0.39200481365977291,1.5660864269415586,0.96661100600682737,1.3138117044789994,0.0968510660401837,-0.41497045683521272,-0.23555261515678594,-1.1177731709718628,-1.5110656097963842,-2.2673569066168415,-2.44782880511669,-2.3097581056178806,-2.1505720078425177,-2.1379733272311197,-3.0522636213910403,-3.0083756662379426,-3.2100022300283029,-3.3367885224363887,-2.4609504973696561,-2.6035414238025734,-1.5896542311152022,-2.23740386452172,-1.3153777123413177,-1.7115592522548193,-0.91303227947332333,-0.19395513239682305,-0.56730752694728914,-0.033725446167341674,-1.1377952255177435,0.24651330440413285,0.86068119701291923,1.0631289708885621,2.0495165788974519,1.5577409978883137,1.5107016521429064,1.8492360562286596,0.96985304639152392,1.3615682390027986,1.7755707232152542,1.2974536303606474,0.062011278800511382,1.4283196952420762,0.15879179903039276,-0.037521341154554633,0.53899592731657164,-0.074088983909330774,-0.62453544488232426,-1.3026221424388009,-1.089978259258485,-2.1191576357264141,-1.0686248423492055,-1.4716336600389677,-2.2368095108575989,-1.790133383775963,-0.28434563396156443,-1.6542774953946853,-1.5300923494764704,-0.838088282844372,-1.3156754366295493,-1.3236271431722795,-1.069706435624554,-0.68623149598402489,-0.018251876740953188,-0.085495541627020388,0.68653418728912141,0.39313701255287842,1.2120952170042669,0.88447522055668193,1.9042426405936714,1.2033440584454427,1.84957484963817,2.0757735161230872,3.4431754286112906,2.20308872891275,2.2093014355255027,3.031900294324871,2.3116252752428434,2.3548012394711706,3.0269199187573745,2.7973553834951783,3.143684942308437,2.4387777818818792,2.705094252970623,2.0956900161959107,2.1950838998829987,1.1139239815440345,-0.33319944167072346,0.044398215722526391,-0.0925196133912504,-0.192125032627519,-0.82057037043538872,-0.59359316003046614,-0.975758711545257,-0.32620388401043532,-2.0977315830039247,-1.5710505609077559,-1.8263517631084056,-1.0349815091462538,-1.3779299325007943,-1.1788569791246857,-1.5331943458098829,-0.88018345684125854,-0.40248642662165224,-1.3139556626533038,-1.2032013149956369,-0.489058666954583,-0.63604468180761131,0.55355986812289781,0.021055643540040664,-0.18982813393673126,0.583173350099703,0.30144881579192773,0.0803876802843837,1.0903672467554253,0.64838662019389748,1.1149595614751862,1.4623491665051644,0.76159233112485836,0.33200754125549276,0.64053819755486519,0.70523061006506371,0.2313957288080829,-0.26586305377421532,-0.56771372293815248,-1.3394038193540618,-1.2033811446415599,-1.4586959327770486,-1.4283596914567067,-1.6121540560358423,-2.1443469709523351,-2.8730266172308645,-3.3392413433297357,-2.7363851844023759,-2.9764716202122146,-2.8308461172506068,-4.0678954496851354,-2.510223529382122,-2.5416387068729067,-2.5671336497294557,-2.49592595214344,-1.7307075118554061],[-1.7452179426654468,-1.8657488130576076,-1.1491298569650248,-1.2813015031003816,-0.84881753482935607,0.6364220531262158,-0.74787827908395677,-0.13178589835027377,2.0705176746146763,1.0309083221603337,1.3884869920572243,2.043421273748991,0.534956610432015,2.4933198622565897,2.2180490250665512,1.2083510551541852,1.6106330852108903,1.9538104745522167,1.2682912271675488,1.3276966229990981,1.3192321472397726,1.424530374068766,1.0861353240694018,1.6211909906712791,1.2791052632182203,1.8379462750452218,0.69302745397147025,1.4879715468464008,1.2204063510808987,1.2338282296790424,0.84375288599512654,0.84118756968960284,0.22404514418561661,1.5023102522125997,0.85817386598007683,1.2000978668511388,1.4325290372130022,0.9986704816072236,1.2730871580806593,1.0382054331929038,1.9020931679484985,0.91460039125115444,2.0604990919928357,1.9557217035439456,2.1794959626250572,1.5709713462672583,2.1131592557271075,1.2975936485014854,1.1599543122997573,0.91211210545588561,1.2113535275874769,0.61124866035062053,0.60104608871000631,0.294860923115794,0.088503364888572925,-1.1085804539311275,-0.41258128213064554,-1.4037990476378375,-0.64371235561709261,-0.68528887183023435,-1.6444452604772597,-1.1755172919196815,-1.5132962963297134,-1.9026094534876208,-1.1839383548806421,-1.7970604489559394,-2.2843442716944451,-1.9665787883694925,-2.0017019381272272,-1.454083665508429,-2.9210607240504149,-1.7527534818684751,-1.8694578645846534,-2.2430151086678456,-2.4086172509721484,-1.7861892445722489,-1.1274495670327,-0.59679886844932906,-0.17381970533451152,-0.37316428037405858,-0.03388513933587721,0.70946987737185219,1.065694677279994,1.5117008395789522,0.19913499782465005,0.75424280011857692,1.6787853150951384,2.2035008797943409,1.3337597806925585,1.6999858485107329,1.6389571410284327,1.81341626422716,1.5597668550391206,2.0644734713976858,0.77298592462464955,2.3870469614150238,0.49267389604806222,1.451274116278559,1.277156516971151,0.351265362400861,0.25080347620247734,-0.34135463409251288,-0.00041758553844956081,-0.5013414534949161,-0.80190427446348844,-0.72540400316086218,-1.4131303976040508,-0.95879767892901346,-1.73015973157817,-1.4980747021779413,-1.8984065098288678,-1.7355442127593999,-1.9229020431196895,-1.7040952804439646,-2.1592215702639117,-1.2690644391971238,-1.9263935083066714,-1.4368328491570743,-2.0292522035345235,-0.59335569409204891,-0.66901343413122949,-1.0410579954407155,-1.1286492469346525,-0.66854622109558393,-0.0871243604235051,0.12961959714512034,0.43384576873884872,-0.013519109159771336,0.43100856006297045,0.45434440818732591,0.33674774088658116,1.1474399925616798,0.97126074569625154,0.66183580494497485,0.65230469994050388,0.83917856825705184,0.49093009431359452,0.36344542266693519,0.73698755091737467,0.52983337439449674,0.53287754367682,0.85444843693919126,1.4640258396516495,0.14265406574038203,1.4603719824261154,0.57950403893108848,0.17716536419444068,0.86068277072203048,0.77116580098300525,0.36254896286223731,0.77377622782177125,0.088963474396281272,-0.3629172469758063,1.1960222830844096,0.64030295869485954,1.6456706511432124,1.7806549603042945,2.139660218487081,2.1009289013155596,2.3023608489764427,1.903213049400529,3.2701498643528071,2.3690844583487407,2.4940753010547514,2.9831292563562242,1.7976820911608247,3.7954613962380543,1.5944746868258535,2.8002726751348814,2.599957304109882,1.4055377911491178,0.52386901262295837,0.78832823868251789,1.0292948710374086,0.85785144501905153,0.43620626937002954,-0.32255827337574361,-0.49452780454188894,0.044596597768497781,-1.0211931031496424,-1.7212269945404362,-1.9740791763704981,-1.8543458439849234,-2.8166977764790651,-2.5347921990042277,-3.6150197469421719,-3.7898952317941625,-4.0126439122492066,-4.65227449536689,-3.1133241268518552,-3.24259640386508,-4.1894910127678457,-3.5192349940805192,-2.9939376340126307,-3.9422719858819892,-3.8054082937073028,-3.2871710797815212,-2.3222887120886337,-2.3696021027271268,-2.2724202129067979],[-2.5411204501867104,-3.606461416340609,-3.5231879201556682,-2.9047193837169609,-2.6181083698008982,-2.5139701870065179,-2.9863505744890748,-2.5408318026094721,-2.26566656310029,-1.8279834924155061,-1.6153072211845885,-1.3368550894161275,-1.5642468980322359,-0.6465665647294685,-0.13370070027009229,-0.48633931230045879,-0.28774150479043636,0.40519077548854943,0.0018312664135661794,1.5387338965659518,0.51888471819947046,1.1538259792051158,0.9655437321138447,1.3299966975935671,1.6549938316758321,2.4022332165684688,1.8227188091342472,2.1261988117410762,2.1535873544980362,2.9029318938304947,1.3436434964864854,2.6548030653774481,2.4129028551936345,1.6255798106875581,2.3295246908557576,2.2202185984936187,2.1677953251388447,2.3984443598582046,1.9923585577991199,2.2125672009554371,2.5985801266452411,1.9498671429771091,1.3978193148182432,0.95529169670198577,2.0775022263027259,0.81123331333547666,1.3717332577208219,1.1210297368198336,0.87974349921848138,0.58851935925717458,-0.034520159343715551,0.49601607763814259,-0.078556731887566345,-0.13700105980230493,-0.28167925064297206,-0.60285656106717367,-1.0752707438139255,-0.78413978798049233,-0.042932370777484685,0.26175653931526621,-0.8305517868632436,-0.55576597012605788,-0.7090400776042256,-1.7976875426034067,-0.91289648460407447,-1.4908607562248535,-1.8501822663098193,-1.4158607617956529,-0.22690902979325012,-1.0156059272306339,-1.5833598231277377,-0.664300408714075,-0.7440012389631141,-1.1033077616249725,-0.61895459551390641,-0.98228397359372188,-1.4105312412530178,-1.3272186917173214,-0.68389385508795175,-1.0048520713348144,-0.69119447096617526,-0.30010064005541748,-0.14662569752228582,-0.89007162023105457,0.16403103010571596,0.077465296222013369,0.56322208452836708,0.78500556926512521,-0.18719558675727632,0.89854955753646837,0.7643914483265557,1.4335084943841814,0.96587782121084254,0.58079695502453965,-0.0061690282206481317,0.56129194404444693,0.50627850766072746,1.2116553377874129,0.46870243910357967,0.43330260045943425,1.176883402410569,2.0143477546158128,0.56509717448993713,0.12081880877269024,-0.33804529368254621,0.26605740672899653,0.85116091856323273,0.351877045416635,0.29483563293156628,0.45802661965230723,0.67753718847956157,-1.1983925109560774,-0.63418028673347571,-0.969937607156156,-1.1299143872493604,-2.2590131218031111,-1.8016986350090463,-0.71667036530356565,-1.1318937519849337,-1.5753996297117412,-1.3951372766757706,-2.5524363532774568,-2.0402399111217959,-2.2249940022709485,-1.7196163158999842,-1.3874603400631718,-2.5766790028536448,-1.5858863718294585,-2.0666444249679206,-1.3991190546352135,-1.4680528712067793,-2.3281664788802767,-0.78545184021095782,-2.3251744571422126,-1.3437330406069188,-0.49048422256507718,-1.8507125272771121,-1.0706767569412627,-0.89309857536331561,-0.40823905677456013,0.019836657967748084,-0.19473989060055863,0.38705810140739372,1.4742883564252103,1.3798050894693756,1.6786111595140545,1.2908715591328366,0.86797193193574362,1.5842991137111573,1.9130609912257877,2.253011577318444,2.2644609858281335,1.7204225069469872,1.4962850330102837,1.5674545547517256,2.0646032181667864,2.6668209533773553,3.2465343743577506,2.4689602407328182,3.0058542708046763,1.9440005900632849,3.2229299300110359,2.574114014125692,3.1243219732361185,2.7740749263684732,2.0710801234053258,1.8318294985513308,1.0757454820283914,2.0702468408726662,1.9381219656206012,1.6877452703383564,2.2032101452843538,1.5887771154767203,1.7771958958879588,1.6404145677688926,1.2127065449086796,0.44230729525696078,0.97160076095904813,0.26344655458128075,-0.46354460359933547,0.527913004728179,1.0390582176447392,-0.81143039815229268,-0.6335603069015544,-0.25516624710765157,-2.1652620818229993,-2.2862860687493485,-1.6672467935298407,-1.232708206971632,-0.26634193994119171,-0.5466624899254835,-1.8243888945203757,-2.1771711269716549,-2.9322534841629997,-3.66737152863446,-2.9350570194965284,-2.8856527507915408,-2.9916840459928919,-2.2540400313213502,-3.044169838702981],[-2.8046825961951054,-3.2431151125110347,-3.1933684052773472,-4.2169758854012258,-2.9855332777238903,-2.8770018201938345,-3.3371466281329583,-2.5431380541012913,-2.225004060548657,-1.8707741041534889,-1.4635237036204929,-0.90922771910405209,-1.1208611100348569,-0.56139587864369367,-0.51023291596470388,0.277453429242338,-0.23673772410957233,0.83839922101046993,-0.077174427257627132,1.1568994706465925,0.79306524798092315,1.1264552856306731,0.87288904663263545,1.0155608985041089,1.236657186482262,1.4790785288512729,1.4250658544290142,1.070154347825552,1.4125049128331479,0.84784500430710985,1.1360373161239803,1.3351474675594381,-0.34198803066279893,-0.44933075435012554,0.6175915618794765,-0.18873191705092007,0.039557172415561975,0.52067345032986179,-0.61797526675898018,-0.79398638183773562,-0.2426280489169714,-0.14489929707198962,0.099214165909012053,0.19351673313464995,-0.64412902529020089,0.11506756762697,0.26150992706613052,-0.61215006204582745,0.086700010321195675,1.317396148221746,0.56510270320506706,1.412152734734345,0.572686827521705,1.7419355046192397,2.1540560180635193,2.3750649984277343,1.646047431433586,1.3990410649570608,2.2707160934460529,1.8086021083861406,1.7750222518262986,2.2396977913132106,2.3828101017693948,0.84467701267252626,0.82181268515723338,0.86911697070340543,0.90339744323254789,0.96646064351238437,0.29919390471243024,-0.48029984843501394,-0.041171250019689021,-0.51653300786505074,-1.13563899561085,-0.981627005363902,-1.3726461557203111,-0.72310165467524434,-1.8839053072860841,-1.2076797589992521,-1.2874250041752617,-2.1419536810434852,-1.4852031698153716,-1.2819036436226383,-0.77586594472108616,-2.0665769327111811,-1.6092360702037549,-0.88222747542012347,-1.8763208676939431,-1.1669941683994447,-1.0049159581432261,0.558796669066713,0.20304735682752761,2.3209940585923405,1.5130753555330152,1.7808859692778,1.5323902443504265,1.6617528322216926,2.393260637234004,2.7014968578677285,2.9214831027109067,2.9122941145141459,2.7551341444680357,2.1319064364213971,2.0189931838638913,1.4346114755693697,1.6065318446370678,1.6310402488403877,0.46642864840693,0.65700636992239814,0.33456715076035171,0.24704312266798545,-0.30994315748124335,-0.65606504916948594,-1.4913142170954408,-1.7461842366743694,-2.1249411278080443,-1.6093523564077885,-1.8027367419279379,-3.1756912019662917,-2.9822072553408105,-2.7444762463716552,-2.3837377158948696,-2.6014866111436219,-2.2817125764506945,-3.3505709525336123,-2.9228195613385717,-3.1520718755775587,-2.2784489213753329,-2.5466605453166804,-1.5485847817369776,-2.544658052196902,-1.4154516220204114,-2.2843218875246123,-1.4026073686988056,-0.735537241509147,-0.25412945246228558,-0.05809133288261889,-0.89257087557513692,0.41072533433272329,0.73665508724299,0.072208922110585716,1.084458592507445,0.37020700021565922,0.27266034020914026,0.80662911272735394,1.0315144639561145,0.90042540845614,0.76390147996869429,0.84024084598368221,0.3068197671714985,0.71583998229141432,0.59244901462561983,-0.586999572222782,-0.17660705719133385,0.51767004330130884,-0.18522955495114068,-0.34141323391116096,-0.16652999825036538,0.44830685680878762,-0.59321969214816261,0.505748735837358,0.62509496307756951,0.61837392242737454,0.20286324045749582,0.15286908470294036,1.3820419694446864,1.3574948444459809,1.4503362868266583,1.7223909836645683,2.0785337816439644,2.294060425657777,2.9339376417250596,2.8724935115553181,2.766208874102916,3.3633126973168874,3.9381939378434847,3.1924309480175364,3.5937593104484185,2.5407645106641472,3.5995697409230107,2.6307935547898555,2.5866859293602928,2.8899770544292225,2.0914691259772269,0.91001233437354179,0.67376904492159606,0.607199323766462,1.1789820261766422,0.24206842673947351,-0.10362579196564112,-0.91037320411928957,0.13236890480897578,-1.3531476313902224,-2.0449451309106341,-1.8985884202709458,-3.1699201043772427,-2.9652250491788315,-2.9332348339165191,-2.1704164449574384,-3.2487845113801606,-3.3720418854559098],[-0.89888439825021293,-0.57619013800552754,-1.146769201116185,-0.46541671432138376,-1.2448492239951185,0.090237280253472968,-0.46612807957583097,0.25540251845255929,0.65837882764014621,0.32740952434974147,1.5535166305824646,0.78747125215251867,0.82014658015736464,1.5237679658829173,1.848159587931641,1.2364315376980506,0.53517249935767608,1.4525299715636559,0.38349714001969848,-0.46050254181896033,-0.18601469791958569,0.81315168262158344,0.782420048507933,-0.44349007266293755,-0.81588124540259954,-2.0623982627274926,-2.3726136337079273,-1.3856725716580032,-2.1489808894101188,-2.6482795255434088,-2.7174542161709394,-2.8501653040323602,-2.7880662361311521,-3.3362173557742287,-2.6983634840177717,-2.625142443912408,-2.2616119039912772,-2.7852284778843797,-1.5519461547297211,-2.2805149838033492,-1.5755962340335832,-1.7683069547614338,-0.35428767521619131,-1.0944093428063586,-0.027237030399535789,-0.2996589261670134,0.71813751792093083,1.147989923924897,0.86140026604722442,1.0883148423683124,1.2571259552987679,1.68882776532828,1.7725699545952793,1.2082249864398715,2.1245095414467023,2.4686416855927851,1.9600181116234294,1.9794864820504863,1.9169518539247647,2.076260788389702,1.5240857153239593,1.1142753414455708,2.284001938845972,0.98812459731350155,1.2778337782904581,0.15426923656270242,1.4620060619301731,0.21969033301477828,0.55342157853688867,0.13372024043443109,-0.45177788194457491,-0.45356793636522358,-1.0486511775635132,-0.7677779335734094,-0.44082556430765207,-1.1376153614646793,-1.1126269808114282,-0.845908989730999,-0.23933434068230086,-0.95186189036723978,-0.79267167560176421,-0.1937796222007003,0.058698786535231889,0.78968804100642831,1.2438489821338492,0.43217590295328878,0.26584528519494588,0.67066829483362733,1.0731268823273967,1.215633628764738,1.9706839900549129,1.9233054397840266,3.0042513881012884,2.5967675324916359,2.192597361916659,1.5152630509462179,2.0936394747639464,2.0218428463173153,1.8602563904056122,1.2029779192495385,0.599120092515724,1.5671056005689046,0.16796209808624096,0.10759313852930202,0.066726797593898712,-0.63406870071519617,-1.3917410574102751,-0.53669125109766547,-3.0727747862547874,-1.7408245526140311,-1.4849085778118933,-2.5264040100660505,-2.692899123137777,-3.07374537418885,-1.6270634906647048,-2.0580630839156209,-2.8601984856415639,-2.8942085277108154,-3.0517033448311852,-1.7040354408355431,-3.0659718154111419,-2.4507690239020552,-1.9629113621492786,-1.9925194928565282,-1.7472924772068907,-0.46508186320664247,-0.53597134063512519,-0.43432465333493231,-0.024767148702625946,0.32278791099158777,0.14640006707387293,0.72613065931529708,0.57439441705304728,1.5400013327754094,1.4672044974023741,1.090494985933754,1.4256156064385979,0.78560497373662164,0.81524876030433879,0.68004903779554093,0.9294604927729746,0.65835647632989813,-0.60031446814701772,-0.31156363704980222,-0.31066246074910669,-0.42244484190452863,-0.81526383878946018,-0.46341450013541918,-0.35615030544698112,-0.68890491308305768,-1.1564516563834986,-0.93007835067507538,-0.961241047642783,-1.31295191145412,-0.58491848995412443,-1.655073318494956,-1.2543919403036989,-1.5109590559132897,-1.4140462196535568,-1.3603338537771912,-0.99141581974098414,-0.964271470351381,-0.16278004640015087,-0.040120935895543386,1.4256661685635392,1.2452315385514203,0.34107229570667708,2.204800792189991,1.6437127966476222,1.3835422986214612,2.6650413854512367,3.1988927929048452,2.8745865541745586,3.0252681091589797,2.4685961056028165,2.7862349295246052,2.7539218627001629,2.5797975903664709,2.9067314021955921,1.9545743792309058,1.9414988834712181,1.6065834028666746,1.7378601908264615,1.8390732172325119,1.7633517582400859,1.5805759935573418,-0.30622897066958282,-0.49901384387420772,0.52458291782700339,-0.54849567087077289,-1.2727701389464192,-0.52105234509241094,-2.1842996040643614,-1.7822427212167169,-2.1165189542746421,-1.1389006083913373,-0.72388218737281174,-1.2895342472657134,-0.82660701269955994,-0.66787320805883676],[1.4627790169348132,1.0866523534084787,2.0079415738890916,1.9142471072366785,1.3623884236347803,1.7903943643303879,2.9356356696908676,1.6975100830043697,2.48756914508054,1.9612243081702476,2.601652236201808,1.7727133213482116,1.5437413117716325,2.1555088957351369,1.4415152995487737,2.3490614629619655,1.7380131028931103,1.066620645065987,0.8046094146047047,-0.29529822213574142,0.41439723609984358,1.2220551261206993,-0.0060520653133026847,-0.95803937904742775,-0.55576578906964857,-0.62310009209656736,-1.8964307682698816,-1.1988636027539248,-0.9521708954035184,-1.7305006957283215,-2.1456670357655807,-2.6246727990697476,-2.1752699688064103,-2.7793323345164818,-2.4748049651283268,-2.5244707343974633,-2.3215736781120051,-2.0232880159584981,-2.4037898442983856,-2.0772723379981164,-2.8728646751656495,-1.4068623965410469,-2.2872436425785834,-1.8101405159286654,-1.4609145332339053,-0.97521988223844258,-1.7033938121696131,-1.175373944744978,-1.3761516273814975,-1.5283138780117091,-0.915498315690142,-0.90565332349361061,-0.77428110435896524,-1.2127324021828043,-0.90738620507782719,-0.44396428452501957,-0.22967281652891386,-0.68252180064106416,0.47566563990045507,-0.9989494768251721,-1.1292074752510415,0.39610729596083849,-0.83921821209179659,-0.4138702266915098,-0.74646573026562635,-0.75476850681125518,-0.45671034957033985,-0.037117469551557161,-0.048656552559596666,0.756141466777217,0.27635064926533681,1.297738932804547,0.44230774182495369,0.27453985362006539,1.4103115843813978,2.136183812354612,1.5739457781821009,2.5196577618256368,1.9852864805486654,1.9850359689636647,1.8618094514801458,3.1772852536085958,2.1324052942084273,2.7295285771738738,2.5861256619799446,2.6403809381882821,2.9022514384400222,3.0971507458805307,3.0882087888583176,2.4556756365221744,2.9193951881451112,2.8420605438736137,2.3162391861309666,2.0431823523117369,1.2943373916232768,1.77155354674505,2.1363006721016617,1.3685507904090013,0.60599990212488908,0.61509989581685542,0.081805275730200189,-0.057634663674602182,-0.60636411210732655,-0.97277920947327978,-0.4731834860062728,-1.473751479208484,-1.456258601698021,-1.5962133197508326,-1.4984809709796068,-2.3814941419143709,-2.0948220323214581,-2.4447318750277121,-2.7506805989803023,-2.5560282403719072,-2.1639671214455731,-2.7685875903006103,-1.502620130290679,-1.2251600788760526,-2.3815895621817909,-2.0217001923674114,-1.434003711688467,-1.2696824280436407,-0.57908196721307481,-1.1062179603776432,-0.0750255025645129,-1.0978880898676084,-0.458694424417614,-0.16904610199669154,-0.66233716701862289,-0.015782292500592533,-0.064489717669186458,0.58931268611754883,0.50410111234237043,1.2074641450797028,0.80718881743381354,0.36957696168024307,-0.6499279883722624,0.21913090086672296,-0.13749169851851917,1.3934850416335092,0.18130857404034523,0.27023426561642094,0.21862400888676947,-0.095297662930364888,-1.008877640041949,0.90328286512220635,-0.0014101050090970513,-0.60509792293039211,-0.33545446616968172,0.57542151968208277,-0.12880056243226373,-0.087157854869078927,-0.30093481408898737,-0.836659771845863,-0.35883075888609028,-0.35708648468403353,0.20042866213480121,0.336570037322124,0.26458451917776254,0.81803437045749827,0.37250765092641769,0.75118354260663689,0.2787858679137224,0.28831776369964873,0.685507976886962,1.2808496437712995,0.98385510068912907,1.0343543911523687,0.50890529942360607,0.47805463767148704,1.0459766393861198,0.90040908202497327,-0.084976942920719334,-0.28786547894967213,0.8403875086751047,0.17806931489955796,-0.036674067142168129,0.08049167611090649,-0.31545148969321057,-0.30501286298111963,-0.31967743781900138,-1.4929085452136597,-0.79386380162256065,-0.91823870672531993,-0.77423588118739572,-0.87627170176739233,-1.6769912704238583,-0.98817855667164711,-1.5236095209826042,-1.1414292511417989,-0.82291925525138365,-0.28125114292851244,-0.38324987060646132,-0.31462153117219832,-0.97034914539750139,0.0042566088871075491,0.75704432308211822,-0.65506024442089672,0.87165151374412642,-0.35052634898117974]],"sep":10} + diff --git a/examples/additive_structure_example.py b/examples/additive_structure_example.py index b8fa1de..6f0c14c 100644 --- a/examples/additive_structure_example.py +++ b/examples/additive_structure_example.py @@ -1,7 +1,7 @@ -import discretesampling.domain.additive_structure as addstruct -from discretesampling.base.algorithms import DiscreteVariableMCMC, DiscreteVariableSMC import numpy as np import pandas as pd +import discretesampling.domain.additive_structure as addstruct +from discretesampling.base.algorithms import DiscreteVariableMCMC, DiscreteVariableSMC def func(x): diff --git a/examples/decision_tree_example.py b/examples/decision_tree_example.py index 97760fc..38201f1 100644 --- a/examples/decision_tree_example.py +++ b/examples/decision_tree_example.py @@ -1,17 +1,11 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Nov 5 14:28:12 2021 +from sklearn import datasets +from sklearn.model_selection import train_test_split -@author: efthi -""" +import numpy as np from discretesampling.domain import decision_tree as dt from discretesampling.base.algorithms import DiscreteVariableMCMC, DiscreteVariableSMC -from sklearn import datasets -from sklearn.model_selection import train_test_split - -import numpy as np data = datasets.load_wine() diff --git a/examples/rjmcmc_example.py b/examples/rjmcmc_example.py new file mode 100644 index 0000000..aecbbc3 --- /dev/null +++ b/examples/rjmcmc_example.py @@ -0,0 +1,170 @@ +import copy +import numpy as np +from scipy.stats import multivariate_normal + +from discretesampling.base.random import RNG +import discretesampling.domain.spectrum as spec +from discretesampling.base.algorithms.rjmcmc import DiscreteVariableRJMCMC +from discretesampling.base.stan_model import StanModel +from discretesampling.base.algorithms.continuous.RandomWalk import RandomWalk + + +stan_model_path = "examples/stanmodels/mixturemodel.stan" + + +# What data should be passed to the stan model given discrete variable x? +def data_function(x): + dim = x.value + data = [["K", dim], ["N", 20], ["y", [3, 4, 5, 4, 5, 3, 4, 5, 6, 3, 4, 5, 15, 16, 15, 17, 16, 18, 17, 14]]] + return data + + +# What params should be evaluated if we've made a discrete move from x to y +# where params is the params vector for the model defined by x +class continuous_proposal(): + def __init__(self): + # grid search parameters + self.grid_size = 1000 + self.min_vals = [0] + self.max_vals = [np.pi] + self.inds = [0] + + # store these values here so that they can be used when eval() is called (for efficiency, so we don't have to + # re-compute many expensive logprob evaluations) + self.current_discrete = None + self.current_continuous = None + self.proposed_discrete = None + self.proposed_continuous = None + self.forward_logprob = None + self.to_remove = None + + def sample(self, x, params, y, rng=RNG()): + dim_x = x.value + dim_y = y.value + + self.current_discrete = x + self.proposed_discrete = y + self.current_continuous = params + + params_temp = copy.deepcopy(params) + if dim_x > 1: + theta = params_temp[0:(dim_x-1)] # k-simplex parameterised as K-1 unconstrained + mu = params_temp[(dim_x-1):(2*dim_x-1)] + sigma = params_temp[(2*dim_x-1):len(params_temp)] + else: + theta = [] # One component, mixing component is undefined + mu = params_temp[0:dim_x] + sigma = params_temp[dim_x:len(params_temp)] + + if (dim_y > dim_x): + # Birth move + # Add new components + n_new_components = dim_y - dim_x + + if dim_x > 0: + theta_new = rng.randomMvNormal(np.zeros(n_new_components), np.identity(n_new_components)) + mu_new = rng.randomMvNormal(np.zeros(n_new_components), np.identity(n_new_components)) + sigma_new = rng.randomMvNormal(np.zeros(n_new_components), np.identity(n_new_components)) + + mvn = multivariate_normal(np.zeros(n_new_components), np.identity(n_new_components)) + forward_logprob = mvn.logpdf(theta_new) + mvn.logpdf(mu_new) + mvn.logpdf(sigma_new) + else: + # initial proposal + theta_new = [] + theta_logprob = 0.0 + if n_new_components > 1: + theta_new = rng.randomMvNormal(np.zeros(n_new_components-1), np.identity(n_new_components-1)) + mvn_theta = multivariate_normal(np.zeros(n_new_components-1), np.identity(n_new_components-1)) + theta_logprob = mvn_theta.logpdf(theta_new) + + mu_new = rng.randomMvNormal(np.zeros(n_new_components), np.identity(n_new_components)) + sigma_new = rng.randomMvNormal(np.zeros(n_new_components), np.identity(n_new_components)) + + mvn = multivariate_normal(np.zeros(n_new_components), np.identity(n_new_components)) + + forward_logprob = theta_logprob + mvn.logpdf(mu_new) + mvn.logpdf(sigma_new) + self.forward_logprob = forward_logprob + + params_temp = np.concatenate((np.array(theta), theta_new, np.array(mu), mu_new, np.array(sigma), sigma_new)) + self.proposed_continuous = params_temp + elif (dim_x > dim_y): + # randomly choose one to remove + to_remove = rng.randomInt(0, dim_x-2) # really need to sort out simplex vs. real param indexing + # birth of component could happen multiple ways (i.e. different n_new_components) + # so I think the reverse_logprob will only be approximate - seems like there might + # be a proof for the summed pdf values of the series of Gaussians that we can use? + mvn = multivariate_normal(0, 1) + if dim_y > 1: + theta = np.delete(theta, to_remove) + else: + theta = np.array([]) + + mu = np.delete(mu, to_remove) + sigma = np.delete(sigma, to_remove) + + forward_logprob = np.log(1/dim_x) + self.forward_logprob = forward_logprob + self.to_remove = to_remove + + params_temp = list(np.concatenate((theta, mu, sigma))) + + self.proposed_continuous = params_temp + return list(params_temp) + + def eval(self, x, params, y, p_params): + dim_x = x.value + dim_y = y.value + + move_logprob = None + + # first check that we're doing an eval for the previous proposal + if self.current_discrete == x: + if self.proposed_discrete == y: + if np.array_equal(self.current_continuous, params): + if np.array_equal(self.proposed_continuous, p_params): + if (dim_y > dim_x): + # Birth move + reverse_logprob = 1 / dim_y + elif (dim_x > dim_y): + # Death move + # Find probabilities of sampling parameters of component that has been removed + n_new_components = dim_x - dim_y + + theta_old = params[self.to_remove] + mu_old = params[dim_x + self.to_remove] + sigma_old = params[2*dim_x + self.to_remove] + + mvn = multivariate_normal(np.zeros(n_new_components), np.identity(n_new_components)) + reverse_logprob = mvn.logpdf(theta_old) + mvn.logpdf(mu_old) + mvn.logpdf(sigma_old) + move_logprob = reverse_logprob - self.forward_logprob + + if move_logprob is None: + raise Exception("Proposal probability undefined. This can only be called directly after sampling a from the \ + proposal.") + + return move_logprob + + +# initialise stan model +model = StanModel(stan_model_path) + +rjmcmc = DiscreteVariableRJMCMC( + spec.SpectrumDimension, + spec.SpectrumDimensionTarget(3, 2), + model, + data_function, + continuous_proposal, + RandomWalk, + False, + False, + True, + 0.5, + spec.SpectrumDimensionInitialProposal(10), +) + +n_chains = 4 +samples = [rjmcmc.sample(10) for c in range(n_chains)] + +dims = [[x[0].value for x in chain] for chain in samples] + +print("Done") diff --git a/examples/sonar_example.py b/examples/sonar_example.py new file mode 100644 index 0000000..dabeca0 --- /dev/null +++ b/examples/sonar_example.py @@ -0,0 +1,143 @@ +import copy +import json +import numpy as np +from discretesampling.base.random import RNG +import discretesampling.domain.spectrum as spec +from discretesampling.base.algorithms.rjmcmc import DiscreteVariableRJMCMC +from discretesampling.base.algorithms.continuous_proposals import sample_offsets, forward_grid_search, reverse_grid_search +from discretesampling.base.stan_model import StanModel +from discretesampling.base.algorithms.continuous.NUTS import NUTS + +stan_model_path = "examples/stanmodels/linear_array.stan" +data_path = "examples/5_targets_noisy.data.json" + +stationary_data = [] + + +# What data should be passed to the stan model given discrete variable x? + + +def data_function(x): + dim = x.value + data = [["M", dim]] + for d in stationary_data: + data.append(d) + return data + + +# What params should be evaluated if we've made a discrete move from x to y +# where params is the params vector for the model defined by x +class continuous_proposal(): + def __init__(self): + # grid search parameters + self.grid_size = 1000 + self.min_vals = [0] + self.max_vals = [np.pi] + self.inds = [0] + + # store these values here so that they can be used when eval() is called (for efficiency, so we don't have to + # re-compute many expensive logprob evaluations) + self.current_discrete = None + self.current_continuous = None + self.proposed_discrete = None + self.proposed_continuous = None + self.forward_logprob = None + self.to_remove = None # don't think there's any way to work out the reverse proposal without this? + + def sample(self, x, params, y, rng=RNG()): + + dim_x = x.value + dim_y = y.value + + self.current_discrete = x + self.proposed_discrete = y + self.current_continuous = params + + params_temp = copy.deepcopy(params) + + if dim_x > 0: + theta = params_temp[0:dim_x] + delta2 = params_temp[dim_x] + else: + theta = np.array([]) + # initialise delta2 + delta2 = rng.randomNormal(0, 1) + params = np.asarray([delta2]) + + if (dim_y > dim_x): + # Birth move + # Add new components + offsets = sample_offsets(self.grid_size, self.min_vals, self.max_vals, rng) + [params_temp, forward_logprob] = forward_grid_search(data_function, model, self.grid_size, self.min_vals, + self.max_vals, offsets, self.inds, params, y) + + elif (dim_x > dim_y): + # randomly choose one to remove + to_remove = rng.randomInt(0, dim_x-1) + if dim_y > 0: + theta = np.delete(theta, to_remove) + else: + theta = np.array([]) + forward_logprob = 1 / dim_x + params_temp = np.hstack((theta, delta2)) + self.to_remove = to_remove + + self.proposed_continuous = params_temp + self.forward_logprob = forward_logprob + + return params_temp + + def eval(self, x, params, y, p_params): + move_logprob = None + dim_x = x.value + dim_y = y.value + # first check that we're doing an eval for the previous proposal + if self.current_discrete == x: + if self.proposed_discrete == y: + if np.array_equal(self.current_continuous, params): + if np.array_equal(self.proposed_continuous, p_params): + if (dim_y > dim_x): + # Birth move + # Add new components + reverse_logprob = 1 / dim_y + elif (dim_x > dim_y): + reverse_logprob = reverse_grid_search(data_function, model, self.grid_size, self.min_vals, + self.max_vals, self.inds, params, x, self.to_remove) + move_logprob = reverse_logprob - self.forward_logprob + + if move_logprob is None: + raise Exception("Proposal probability undefined. This can only be called directly after sampling a from the \ + proposal.") + + return move_logprob + + +# initialise stan model +model = StanModel(stan_model_path) + +rjmcmc = DiscreteVariableRJMCMC( + spec.SpectrumDimension, + spec.SpectrumDimensionTarget(3, 2), + model, + data_function, + continuous_proposal, + NUTS, + False, + False, + True, + 0.5, + spec.SpectrumDimensionInitialProposal(1), + 5 +) + +with open(data_path, 'r') as f: + json_data = (json.load(f)) + for key, value in json_data.items(): + stationary_data.append([key, value]) + +n_chains = 1 +samples = [rjmcmc.sample(10) for c in range(n_chains)] + +dims = [[x[0].value for x in chain] for chain in samples] + +print("Done") diff --git a/examples/sonar_example_smc.py b/examples/sonar_example_smc.py new file mode 100644 index 0000000..b9502a0 --- /dev/null +++ b/examples/sonar_example_smc.py @@ -0,0 +1,224 @@ +import copy +import json +import numpy as np +from discretesampling.base.random import RNG +from discretesampling.base.algorithms.continuous.RandomWalk import RandomWalk +import discretesampling.base.reversible_jump as rj +import discretesampling.domain.spectrum as spec +from discretesampling.base.algorithms import DiscreteVariableSMC +from discretesampling.base.algorithms.continuous_proposals import ( + sample_offsets, grid_search_logprobs, forward_grid_search, reverse_grid_search +) +from discretesampling.base.stan_model import StanModel + +stan_model_path = "examples/stanmodels/linear_array.stan" +data_path = "examples/5_targets_noisy.data.json" + +stationary_data = [] + + +# What data should be passed to the stan model given discrete variable x? +def data_function(x): + dim = x.value + data = [["M", dim]] + for d in stationary_data: + data.append(d) + return data + + +# What params should be evaluated if we've made a discrete move from x to y +# where params is the params vector for the model defined by x +class continuous_proposal: + def __init__(self): + # grid search parameters + self.grid_size = 100 + self.min_vals = [0] + self.max_vals = [np.pi] + self.inds = [0] + + def sample(self, x, params, y, rng=RNG()): + + dim_x = 0 + dim_y = 0 + if hasattr(x, 'value'): + dim_x = x.value + else: + dim_x = len(x) + + if y is not None: + if hasattr(y, 'value'): + dim_y = y.value + else: + dim_y = len(y) + + params_temp = copy.deepcopy(params) + + logprob_vals = np.asarray([]) + offset_vals = np.asarray([]) + + if dim_x > 0: + theta = params_temp[0:dim_x] + delta2 = params_temp[dim_x] + else: + theta = np.array([]) + # initialise delta2 + delta2 = rng.randomNormal(0, 1) + params = np.asarray([delta2]) + + if (dim_y > dim_x): + # Birth move + # Add new components + offsets = sample_offsets(self.grid_size, self.min_vals, self.max_vals, rng) + [params_temp, forward_logprob] = forward_grid_search(data_function, model, self.grid_size, self.min_vals, + self.max_vals, offsets, self.inds, params[0:(dim_x + 1)], y) + if len(params) > dim_x+1: + if params[dim_x + 1] == 1: + logprob_vals = params_temp[-(dim_x + 1):-1] + offset_vals = params_temp[(dim_x + 2):-(dim_x + 1)] + elif params[dim_x + 1] == -1: + logprob_vals = params_temp[-(dim_x + 2):-2] # remove one logprob value as it's death -> birth + offset_vals = params_temp[(dim_x + 2):-(dim_x + 3)] + else: + raise Exception("Something weird happened with the discrete moves.") + last_move = 1 # flag birth + # store previous parameters for eval check in case we go through a series of discrete moves that lead back to the + # same point (expensive to work out forward_logprob again) + params_temp = np.hstack((params_temp, last_move, offset_vals, offsets, logprob_vals, forward_logprob)) + + elif (dim_x > dim_y): + # randomly choose one to remove + to_remove = rng.randomInt(0, dim_x-1) + if dim_y > 0: + theta = np.delete(theta, to_remove) + else: + theta = np.array([]) + forward_logprob = 1 / dim_x + if len(params) > dim_x+1: + if params[dim_x + 1] == 1: + logprob_vals = params_temp[-(dim_x + 2):-2] # remove one logprob value as it's birth -> death + offset_vals = params_temp[(dim_x + 2):-(dim_x + 3)] + elif params[dim_x + 1] == -1: + logprob_vals = params_temp[-(dim_x + 2):-3] # remove two logprob value as it's death -> death + offset_vals = params_temp[(dim_x + 2):-(dim_x + 3)] + else: + raise Exception("Something weird happened with the discrete moves.") + last_move = -1 # flag death + params_temp = np.hstack((theta, delta2, last_move, offset_vals, logprob_vals, forward_logprob)) + + return params_temp + + def eval(self, x, params, y, p_params): + dim_x = x.value + dim_y = y.value + current_data = data_function(x) + param_length = model.num_unconstrained_parameters(current_data) + p_data = data_function(y) + p_param_length = model.num_unconstrained_parameters(p_data) + if dim_x == 0 and dim_y == 1: + # initialisation (birth move) + reverse_logprob = 1 / dim_y + move_logprob = reverse_logprob - p_params[-1] + elif dim_x == 0: + # undefined move + move_logprob = np.NINF + elif np.array_equal(p_params[p_param_length:(p_param_length + param_length)], params[0:param_length]): + if (dim_y == dim_x + 1): + # Birth move + num_matching = 0 + for n in range(0, param_length): + num_matching += len(np.where(p_params == params[n])) + if num_matching == param_length: + # can reuse forward logprob + reverse_logprob = 1 / dim_y + move_logprob = reverse_logprob - p_params[-1] + else: + # requires multiple moves + move_logprob = np.NINF + elif (dim_y == dim_x + 1): + # Death move + num_matching = 0 + for n in range(0, p_param_length): + num_matching += len(np.where(params == p_params[n])) + if num_matching == p_param_length: + # can reuse forward logprob + all_inds = range(0, p_param_length) + for n in range(0, p_param_length): + if len(np.where(params == all_inds[n])) == 0: + death_ind = all_inds[n] + reverse_logprob = reverse_grid_search(data_function, model, self.grid_size, self.min_vals, + self.max_vals, self.inds, params, x, death_ind) + move_logprob = reverse_logprob - p_params[-1] + else: + # requires multiple moves + move_logprob = np.NINF + else: + # requires multiple moves + move_logprob = np.NINF + + return move_logprob + + def eval_all(self, x, params, y, p_params): + # evaluate all possible + dim_x = x.value + dim_y = y.value + + evals = [] + if (dim_y == dim_x + 1): + # Birth move + reverse_logprob = 1 / dim_y + # calculate move probability for grid specified by offsets sampled in p_params + offsets = np.asarray([p_params[-dim_y-1]]) + [proposals, forward_logprobs] = grid_search_logprobs(data_function, model, self.grid_size, self.min_vals, + self.max_vals, offsets, self.inds, params[0:(dim_x + 1)], y) + for proposal, forward_logprob in zip(proposals, forward_logprobs): + move_logprob = reverse_logprob - forward_logprob + evals.append((proposal, move_logprob)) + elif (dim_y == dim_x - 1): + # Death move + theta0 = params[0:dim_x] + delta2 = params[dim_x] + forward_logprob = 1 / dim_y + for death_ind in range(0, dim_x): + if dim_y > 0: + theta = np.delete(theta0, death_ind) + else: + theta = np.array([]) + reverse_logprob = reverse_grid_search(data_function, model, self.grid_size, self.min_vals, + self.max_vals, self.inds, params, x, death_ind) + move_logprob = reverse_logprob - forward_logprob + evals.append((np.hstack((theta, delta2, params[dim_x + 1:-1])), move_logprob)) + else: + evals = [(np.asarray([]), np.NINF)] + + return evals + + +# initialise stan model +model = StanModel(stan_model_path) + +# set variables used in the proposal +rj.set_proposal_attributes( + spec.SpectrumDimensionTarget(3, 2), + model, + data_function, + continuous_proposal, + RandomWalk, + 0.5 +) + +# define target +target = rj.ReversibleJumpTarget() +initialProposal = rj.ReversibleJumpInitialProposal(spec.SpectrumDimension, spec.SpectrumDimensionInitialProposal(1)) + +specSMC = DiscreteVariableSMC(rj.ReversibleJumpVariable, target, initialProposal) + +with open(data_path, 'r') as f: + json_data = (json.load(f)) + for key, value in json_data.items(): + stationary_data.append([key, value]) + +try: + SMCSamples = specSMC.sample(5, 10) + +except ZeroDivisionError: + print("SMC sampling failed due to division by zero") diff --git a/examples/stan_model_example.py b/examples/stan_model_example.py new file mode 100644 index 0000000..e82ad15 --- /dev/null +++ b/examples/stan_model_example.py @@ -0,0 +1,43 @@ +import math +import numpy as np +from scipy.stats import multivariate_normal +from discretesampling.base.random import RNG +from discretesampling.base.stan_model import StanModel + + +mixturemodel = StanModel( + "examples/stanmodels/mixturemodel.stan" +) + +niters = 500 +sigma = 1 + +data = [["K", 5], ["N", 20], ["y", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]]] + +# 5 mixture component locations and sigmas, 5-simplex (described by four-element unconstrained vector) +param_length = mixturemodel.num_unconstrained_parameters(data) +print("Param_length:", param_length) + +mu = [0 for i in range(param_length)] +sigma = np.identity(param_length) + +rng = RNG(seed=0) +init = rng.randomMvNormal(mu, sigma*10) +samples = [init] +current = init + +# MCMC +for i in range(niters): + proposed = current + rng.randomMvNormal(mu, sigma) + current_target = mixturemodel.eval(data, current)[0] + proposed_target = mixturemodel.eval(data, proposed)[0] + forward_logprob = multivariate_normal.logpdf(proposed, mean=current, cov=sigma) + reverse_logprob = multivariate_normal.logpdf(current, mean=proposed, cov=sigma) + log_acceptance_ratio = proposed_target - current_target + reverse_logprob - forward_logprob + if log_acceptance_ratio > 0: + log_acceptance_ratio = 0 + acceptance_probability = min(1, math.exp(log_acceptance_ratio)) + q = rng.random() + if (q < acceptance_probability): + current = proposed + samples.append(current) diff --git a/Stan for RJMCMC problems/decayingsinusoids.stan b/examples/stanmodels/decayingsinusoids.stan similarity index 100% rename from Stan for RJMCMC problems/decayingsinusoids.stan rename to examples/stanmodels/decayingsinusoids.stan diff --git a/examples/stanmodels/linear_array.stan b/examples/stanmodels/linear_array.stan new file mode 100644 index 0000000..fd87fa7 --- /dev/null +++ b/examples/stanmodels/linear_array.stan @@ -0,0 +1,82 @@ +data { + int M; // number of signals + int L; // number of sensors + int T; // number of samples for each sensor + real fs; // sampling frequency + real wave_speed; // wave speed + array[L] vector[T] d; // samples from the sensors + real sep; // separation between sensors +} + +transformed data{ + int N = T*L; + int f_n = T/2; // Nyquist frequency + real N_2 = N/2.0; + real dTd; + array[f_n, L, 2] real GHd_factor; + real sinfunc; + array[f_n] real sumsinsq; + array[f_n] real sumcossq; + array[f_n] real sumsincos; + + dTd = 0.0; + for (l in 1:L){ + dTd = dTd + sum(d[l].*d[l]); + } + for (l in 1:L){ + GHd_factor[:, l, 1] = rep_array(0.0, f_n); + GHd_factor[:, l, 2] = rep_array(0.0, f_n); + } + for (omega in 1:f_n){ + // calculate Fourier transform + for (t in 0:(T-1)){ + for (l in 1:L){ + GHd_factor[omega, l, 1] = GHd_factor[omega, l, 1] + sin(omega*t/fs)*d[l,t+1]; + GHd_factor[omega, l, 2] = GHd_factor[omega, l, 2] + cos(omega*t/fs)*d[l,t+1]; + } + } + // calculate other dot products + sinfunc = sin(4*pi()*omega/fs)/(8*pi()*omega/(T*fs)); + sumsinsq[omega] = T/2 - sinfunc; + sumcossq[omega] = T/2 + sinfunc; + sumsincos[omega] = (1-cos(4*pi()*omega/fs))/(8*pi()*omega/(T*fs)); + } +} + +parameters { + vector[M] theta; // directions of arrival for each signal + real delta2; // signal to noise ratio +} + +model { + matrix[M,M] GHG; + vector[M] GHd; + real delta2_const; + real alpha_factor; + array[M, L] real alpha; + delta2 ~ inv_gamma(2,0.1); + delta2_const = 1+1/(N*delta2); + for (omega in 1:f_n){ + GHd = rep_vector(0.0, M); + GHG = rep_matrix(rep_row_vector(0.0, M), M); + for (m in 1:M){ + GHG[m,m] = N; + alpha_factor = 2*pi()*(omega/fs)*sep*cos(theta[m])/(wave_speed); + for (l in 1:L){ + alpha[m,l] = (l-1)*alpha_factor; + GHd[m] = GHd[m] + cos(alpha[m,l])*GHd_factor[omega,l,1] + sin(alpha[m,l])*GHd_factor[omega,l,2]; + } + } + for (m in 1:M){ + for (m1 in m+1:M){ + for (l in 1:L){ + GHG[m,m1] = GHG[m,m1] + sin(alpha[m,l])*sin(alpha[m1,l])*sumcossq[omega] + (cos(alpha[m,1])*sin(alpha[m1,l])+cos(alpha[m1,l])*sin(alpha[m,l]))*sumsincos[omega] + cos(alpha[m,l])*cos(alpha[m1,l])*sumsinsq[omega]; + } + GHG[m,m1] = delta2_const*GHG[m,m1]; + GHG[m1,m] = GHG[m,m1]; + } + } + target += -N_2*log((dTd - (GHd'/GHG)*GHd)/2); + } + target += -f_n*log(N*delta2 + delta2_const)/2 +2200000; +} \ No newline at end of file diff --git a/Stan for RJMCMC problems/mixturemodel.data.json b/examples/stanmodels/mixturemodel.data.json similarity index 100% rename from Stan for RJMCMC problems/mixturemodel.data.json rename to examples/stanmodels/mixturemodel.data.json diff --git a/Stan for RJMCMC problems/mixturemodel.stan b/examples/stanmodels/mixturemodel.stan similarity index 100% rename from Stan for RJMCMC problems/mixturemodel.stan rename to examples/stanmodels/mixturemodel.stan diff --git a/Stan for RJMCMC problems/probabilistic_tree.stan b/examples/stanmodels/probabilistic_tree.stan similarity index 100% rename from Stan for RJMCMC problems/probabilistic_tree.stan rename to examples/stanmodels/probabilistic_tree.stan diff --git a/Stan for RJMCMC problems/tree_structure.r b/examples/stanmodels/tree_structure.r similarity index 100% rename from Stan for RJMCMC problems/tree_structure.r rename to examples/stanmodels/tree_structure.r diff --git a/requirements.txt b/requirements.txt index 6b9ea81..69a2cc9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,10 +3,12 @@ numpy scipy scikit-learn pandas +cmdstanpy +bridgestan mpi4py tqdm flake8 pytest pytest-mpi pytest-cov -coverage \ No newline at end of file +coverage diff --git a/setup.cfg b/setup.cfg index f4a6e02..d7666a5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -7,6 +7,10 @@ url = https://github.com/DiscreteVariablesTaskForce/DiscreteSamplingFramework [flake8] max-line-length = 127 +exclude = + .git + .github + examples/stanmodels [tool:pytest] addopts = diff --git a/setup.py b/setup.py index b024da8..3e1b096 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,6 @@ -from setuptools import setup +from setuptools import setup, find_packages -setup() +setup( + packages=find_packages(exclude=("examples", "tests", "bridgestan",)) +)