From 7953ecb3723c78a5e139041119991cf42d0bcc51 Mon Sep 17 00:00:00 2001 From: aryandeshwal Date: Sun, 6 Jun 2021 19:54:42 -0700 Subject: [PATCH] Upload HyBO code --- GPmodel/__init__.py | 0 GPmodel/inference/__init__.py | 0 GPmodel/inference/inference.py | 114 ++++++ GPmodel/kernels/__init__.py | 0 GPmodel/kernels/mixeddiffusionkernel.py | 102 ++++++ GPmodel/kernels/mixedkernel.py | 33 ++ GPmodel/likelihoods/__init__.py | 0 GPmodel/likelihoods/gaussian.py | 31 ++ GPmodel/likelihoods/likelihood.py | 7 + GPmodel/means/__init__.py | 0 GPmodel/means/constant.py | 31 ++ GPmodel/means/mean.py | 7 + GPmodel/models/__init__.py | 0 GPmodel/models/gp.py | 34 ++ GPmodel/models/gp_regression.py | 21 ++ GPmodel/modules/__init__.py | 0 GPmodel/modules/gp_modules.py | 19 + GPmodel/sampler/__init__.py | 0 GPmodel/sampler/priors.py | 113 ++++++ GPmodel/sampler/sample_edgeweight.py | 62 ++++ GPmodel/sampler/sample_hyper.py | 118 ++++++ GPmodel/sampler/sample_mixed_posterior.py | 102 ++++++ GPmodel/sampler/sample_posterior.py | 88 +++++ GPmodel/sampler/slice_lengthscale.py | 54 +++ GPmodel/sampler/slice_log_order_variance.py | 52 +++ GPmodel/sampler/tool_partition.py | 123 +++++++ GPmodel/sampler/tool_slice_sampling.py | 93 +++++ acquisition/__init__.py | 0 acquisition/acquisition_functions.py | 18 + acquisition/acquisition_marginalization.py | 77 ++++ acquisition/acquisition_optimization.py | 114 ++++++ .../acquisition_optimizers/__init__.py | 0 .../continuous_optimizer.py | 51 +++ .../acquisition_optimizers/graph_utils.py | 50 +++ .../acquisition_optimizers/greedy_ascent.py | 43 +++ .../simulated_annealing.py | 79 ++++ .../acquisition_optimizers/starting_points.py | 62 ++++ config.py | 8 + experiments/__init__.py | 0 experiments/exp_utils.py | 25 ++ experiments/random_seed_config.py | 22 ++ experiments/test_functions/__init__.py | 0 experiments/test_functions/em_func.py | 100 ++++++ .../experiment_configuration.py | 155 ++++++++ experiments/test_functions/mixed_integer.py | 62 ++++ experiments/test_functions/ml_datasets.py | 87 +++++ experiments/test_functions/nn_ml_datasets.py | 126 +++++++ .../test_functions/pressure_vessel_design.py | 93 +++++ experiments/test_functions/push_robot_14d.py | 108 ++++++ experiments/test_functions/push_world.py | 275 ++++++++++++++ .../test_functions/robot_push_14d/__init__.py | 0 .../__pycache__/__init__.cpython-36.pyc | Bin 0 -> 233 bytes .../__pycache__/__init__.cpython-37.pyc | Bin 0 -> 169 bytes .../__pycache__/push_function.cpython-36.pyc | Bin 0 -> 2774 bytes .../__pycache__/push_function.cpython-37.pyc | Bin 0 -> 2606 bytes .../__pycache__/push_utils.cpython-36.pyc | Bin 0 -> 8096 bytes .../__pycache__/push_utils.cpython-37.pyc | Bin 0 -> 7926 bytes .../robot_push_14d/push_function.py | 75 ++++ .../robot_push_14d/push_utils.py | 234 ++++++++++++ .../robot_push_14d/rover_function.py | 253 +++++++++++++ .../robot_push_14d/rover_utils.py | 338 ++++++++++++++++++ .../robot_push_14d/simple_functions.py | 120 +++++++ experiments/test_functions/speed_reducer.py | 121 +++++++ experiments/test_functions/weld_design.py | 126 +++++++ main.py | 188 ++++++++++ requirement.txt | 12 + utils.py | 64 ++++ 67 files changed, 4290 insertions(+) create mode 100755 GPmodel/__init__.py create mode 100755 GPmodel/inference/__init__.py create mode 100755 GPmodel/inference/inference.py create mode 100755 GPmodel/kernels/__init__.py create mode 100644 GPmodel/kernels/mixeddiffusionkernel.py create mode 100644 GPmodel/kernels/mixedkernel.py create mode 100755 GPmodel/likelihoods/__init__.py create mode 100755 GPmodel/likelihoods/gaussian.py create mode 100755 GPmodel/likelihoods/likelihood.py create mode 100755 GPmodel/means/__init__.py create mode 100755 GPmodel/means/constant.py create mode 100755 GPmodel/means/mean.py create mode 100755 GPmodel/models/__init__.py create mode 100755 GPmodel/models/gp.py create mode 100755 GPmodel/models/gp_regression.py create mode 100755 GPmodel/modules/__init__.py create mode 100755 GPmodel/modules/gp_modules.py create mode 100755 GPmodel/sampler/__init__.py create mode 100755 GPmodel/sampler/priors.py create mode 100755 GPmodel/sampler/sample_edgeweight.py create mode 100755 GPmodel/sampler/sample_hyper.py create mode 100644 GPmodel/sampler/sample_mixed_posterior.py create mode 100755 GPmodel/sampler/sample_posterior.py create mode 100644 GPmodel/sampler/slice_lengthscale.py create mode 100644 GPmodel/sampler/slice_log_order_variance.py create mode 100755 GPmodel/sampler/tool_partition.py create mode 100755 GPmodel/sampler/tool_slice_sampling.py create mode 100755 acquisition/__init__.py create mode 100755 acquisition/acquisition_functions.py create mode 100755 acquisition/acquisition_marginalization.py create mode 100755 acquisition/acquisition_optimization.py create mode 100755 acquisition/acquisition_optimizers/__init__.py create mode 100644 acquisition/acquisition_optimizers/continuous_optimizer.py create mode 100755 acquisition/acquisition_optimizers/graph_utils.py create mode 100755 acquisition/acquisition_optimizers/greedy_ascent.py create mode 100755 acquisition/acquisition_optimizers/simulated_annealing.py create mode 100755 acquisition/acquisition_optimizers/starting_points.py create mode 100755 config.py create mode 100755 experiments/__init__.py create mode 100755 experiments/exp_utils.py create mode 100755 experiments/random_seed_config.py create mode 100755 experiments/test_functions/__init__.py create mode 100644 experiments/test_functions/em_func.py create mode 100755 experiments/test_functions/experiment_configuration.py create mode 100644 experiments/test_functions/mixed_integer.py create mode 100644 experiments/test_functions/ml_datasets.py create mode 100644 experiments/test_functions/nn_ml_datasets.py create mode 100644 experiments/test_functions/pressure_vessel_design.py create mode 100644 experiments/test_functions/push_robot_14d.py create mode 100644 experiments/test_functions/push_world.py create mode 100755 experiments/test_functions/robot_push_14d/__init__.py create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/__init__.cpython-36.pyc create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/__init__.cpython-37.pyc create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-36.pyc create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-37.pyc create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/push_utils.cpython-36.pyc create mode 100644 experiments/test_functions/robot_push_14d/__pycache__/push_utils.cpython-37.pyc create mode 100755 experiments/test_functions/robot_push_14d/push_function.py create mode 100755 experiments/test_functions/robot_push_14d/push_utils.py create mode 100755 experiments/test_functions/robot_push_14d/rover_function.py create mode 100755 experiments/test_functions/robot_push_14d/rover_utils.py create mode 100755 experiments/test_functions/robot_push_14d/simple_functions.py create mode 100644 experiments/test_functions/speed_reducer.py create mode 100644 experiments/test_functions/weld_design.py create mode 100644 main.py create mode 100644 requirement.txt create mode 100755 utils.py diff --git a/GPmodel/__init__.py b/GPmodel/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/inference/__init__.py b/GPmodel/inference/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/inference/inference.py b/GPmodel/inference/inference.py new file mode 100755 index 0000000..0544b9d --- /dev/null +++ b/GPmodel/inference/inference.py @@ -0,0 +1,114 @@ +import numpy as np + +import torch +import torch.nn as nn + + +class Inference(nn.Module): + + def __init__(self, train_data, model): + super(Inference, self).__init__() + self.model = model + self.train_x = train_data[0] + self.train_y = train_data[1] + self.output_min = torch.min(self.train_y) + self.output_max = torch.max(self.train_y) + self.mean_vec = None + self.gram_mat = None + # cholesky is lower triangular matrix + self.cholesky = None + self.jitter = 0 + + def gram_mat_update(self, hyper=None): + if hyper is not None: + self.model.vec_to_param(hyper) + + self.mean_vec = self.train_y - self.model.mean(self.train_x.float()) + self.gram_mat = self.model.kernel(self.train_x) + torch.diag(self.model.likelihood(self.train_x.float())) + + def cholesky_update(self, hyper): + self.gram_mat_update(hyper) + eye_mat = torch.diag(self.gram_mat.new_ones(self.gram_mat.size(0))) + for jitter_const in [0, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3]: + chol_jitter = torch.trace(self.gram_mat).item() * jitter_const + try: + # cholesky is lower triangular matrix + self.cholesky = torch.cholesky(self.gram_mat + eye_mat * chol_jitter, upper=False) + self.jitter = chol_jitter + return + except RuntimeError: + pass + raise RuntimeError('Absolute entry values of Gram matrix are between %.4E~%.4E with trace %.4E' % + (torch.min(torch.abs(self.gram_mat)).item(), torch.max(torch.abs(self.gram_mat)).item(), + torch.trace(self.gram_mat).item())) + + def predict(self, pred_x, hyper=None, verbose=False, compute_grad = False): + if hyper is not None: + param_original = self.model.param_to_vec() + self.cholesky_update(hyper) + + k_pred_train = self.model.kernel(pred_x, self.train_x) + k_pred = self.model.kernel(pred_x, diagonal=True) + + # cholesky is lower triangular matrix + chol_solver = torch.triangular_solve(torch.cat([k_pred_train.t(), self.mean_vec], 1), self.cholesky, upper=False)[0] + chol_solve_k = chol_solver[:, :-1] + chol_solve_y = chol_solver[:, -1:] + + pred_mean = torch.mm(chol_solve_k.t(), chol_solve_y) + self.model.mean(pred_x) + pred_quad = (chol_solve_k ** 2).sum(0).view(-1, 1) + pred_var = k_pred - pred_quad + + if verbose: + numerically_stable = (pred_var >= 0).all() + zero_pred_var = (pred_var <= 0).all() + + if hyper is not None: + self.cholesky_update(param_original) + + if compute_grad: + alpha = torch.cholesky_solve(self.mean_vec, self.cholesky, upper=False) + grad_cross = self.model.kernel.grad(self.train_x, pred_x) + grad_xp_m = torch.mm(grad_cross, k_pred_train.t()*alpha) + gamma = torch.triangular_solve(chol_solve_k, self.cholesky.t(), upper=True)[0] + grad_xp_v = -2 * torch.mm(gamma.t(), (grad_cross * k_pred_train).t()).t() + return pred_mean, pred_var.clamp(min=1e-8), grad_xp_m, grad_xp_v + else: + if verbose: + return pred_mean, pred_var.clamp(min=1e-8), numerically_stable, zero_pred_var + else: + return pred_mean, pred_var.clamp(min=1e-8) + + + + def negative_log_likelihood(self, hyper=None): + if hyper is not None: + param_original = self.model.param_to_vec() + self.cholesky_update(hyper) + + # cholesky is lower triangular matrix + mean_vec_sol = torch.triangular_solve(self.mean_vec, self.cholesky, upper=False)[0] + nll = 0.5 * torch.sum(mean_vec_sol ** 2) + torch.sum(torch.log(torch.diag(self.cholesky))) + 0.5 * self.train_y.size(0) * np.log(2 * np.pi) + if hyper is not None: + self.cholesky_update(param_original) + return nll + + +if __name__ == '__main__': + n_size_ = 50 + jitter_const_ = 0 + for _ in range(10): + A_ = torch.randn(n_size_, n_size_ - 2) + A_ = A_.matmul(A_.t()) * 0 + 1e-6 + A_ = A_ + torch.diag(torch.ones(n_size_)) * jitter_const_ * torch.trace(A_).item() + b_ = torch.randn(n_size_, 3) + L_ = torch.cholesky(A_, upper=False) + assert (torch.diag(L_) > 0).all() + abs_min = torch.min(torch.abs(A_)).item() + abs_max = torch.max(torch.abs(A_)).item() + trace = torch.trace(A_).item() + print(' %.4E~%.4E %.4E' % (abs_min, abs_max, trace)) + print(' jitter:%.4E' % (trace * jitter_const_)) + print('The smallest eigen value : %.4E\n' % torch.min(torch.diag(L_)).item()) + torch.triangular_solve(b_, L_, upper=False) + diff --git a/GPmodel/kernels/__init__.py b/GPmodel/kernels/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/kernels/mixeddiffusionkernel.py b/GPmodel/kernels/mixeddiffusionkernel.py new file mode 100644 index 0000000..a6c72f6 --- /dev/null +++ b/GPmodel/kernels/mixeddiffusionkernel.py @@ -0,0 +1,102 @@ +import math + +import torch +from GPmodel.kernels.mixedkernel import MixedKernel +import numpy as np + +class MixedDiffusionKernel(MixedKernel): + def __init__(self, log_order_variances, grouped_log_beta, fourier_freq_list, fourier_basis_list, lengthscales, num_discrete, num_continuous): + super(MixedDiffusionKernel, self).__init__(log_order_variances, grouped_log_beta, fourier_freq_list, fourier_basis_list, lengthscales, num_discrete, num_continuous) + + def forward(self, x1, x2=None, diagonal=False): + """ + :param x1, x2: each row is a vector with vertex numbers starting from 0 for each + :return: + """ + if diagonal: + assert x2 is None + stabilizer = 0 + if x2 is None: + x2 = x1 + if diagonal: + stabilizer = 1e-6 * x1.new_ones(x1.size(0), 1, dtype=torch.float32) + else: + stabilizer = torch.diag(1e-6 * x1.new_ones(x1.size(0), dtype=torch.float32)) + + base_kernels = [] + for i in range(len(self.fourier_freq_list)): + beta = torch.exp(self.grouped_log_beta[i]) + fourier_freq = self.fourier_freq_list[i] + fourier_basis = self.fourier_basis_list[i] + cat_i = fourier_freq.size(0) + discrete_kernel = ((1-torch.exp(-beta*cat_i))/(1+(cat_i-1)*torch.exp(-beta*cat_i)))**((x1[:, i].unsqueeze(1)[:, np.newaxis] != x2[:, i].unsqueeze(1)).sum(axis=-1)) + if diagonal: + base_kernels.append(torch.diagonal(discrete_kernel).unsqueeze(1)) + else: + base_kernels.append(discrete_kernel) + + lengthscales = torch.exp(self.lengthscales)**2 + temp_x_1 = x1[:, self.num_discrete:]/lengthscales + temp_x_2 = x2[:, self.num_discrete:]/lengthscales + + for i in range(self.num_continuous): + normalized_dists = torch.cdist(temp_x_1[:, i].unsqueeze(1), temp_x_2[:, i].unsqueeze(1)) + gaussian_kernel = torch.exp(-0.5 * (normalized_dists) ** 2) + if not diagonal: + base_kernels.append(gaussian_kernel) + else: + base_kernels.append(torch.diagonal(gaussian_kernel).unsqueeze(1)) + base_kernels = torch.stack(base_kernels) + if diagonal: + base_kernels = base_kernels.squeeze(-1) + + num_dimensions = self.num_discrete + self.num_continuous + if (not diagonal): + e_n = torch.empty([num_dimensions + 1, \ + base_kernels.size(1), base_kernels.size(2)]) + e_n[0, :, :] = 1.0 + interaction_orders = torch.arange(1, num_dimensions+1).reshape([-1, 1, 1, 1]).float() + kernel_dim = -3 + shape = [1 for _ in range(3)] + else: + e_n = torch.empty([num_dimensions + 1, \ + base_kernels.size(1)]) + e_n[0, :] = 1.0 + interaction_orders = torch.arange(1, num_dimensions+1).reshape([-1, 1, 1]).float() + kernel_dim = -2 + shape = [1 for _ in range(2)] + + s_k = base_kernels.unsqueeze(kernel_dim - 1).pow(interaction_orders).sum(dim=kernel_dim) + m1 = torch.tensor([-1.0]) + shape[kernel_dim] = -1 + + + for deg in range(1, num_dimensions + 1): # deg goes from 1 to R (it's 1-indexed!) + ks = torch.arange(1, deg + 1, dtype=torch.float).reshape(*shape) # use for pow + kslong = torch.arange(1, deg + 1, dtype=torch.long) # use for indexing + # note that s_k is 0-indexed, so we must subtract 1 from kslong + sum_ = ( + m1.pow(ks - 1) * e_n.index_select(kernel_dim, deg - kslong) * s_k.index_select(kernel_dim, kslong - 1) + ).sum(dim=kernel_dim) / deg + if kernel_dim == -3: + e_n[deg, :, :] = sum_ + else: + e_n[deg, :] = sum_ + + order_variances = torch.exp(self.log_order_variances) + if kernel_dim == -3: + kernel_mat = torch.exp(self.log_amp) * ((order_variances.unsqueeze(-1).unsqueeze(-1) * e_n.narrow(kernel_dim, 1, num_dimensions)).sum(dim=kernel_dim)) + stabilizer + return torch.exp(self.log_amp) * ((order_variances.unsqueeze(-1).unsqueeze(-1) * e_n.narrow(kernel_dim, 1, num_dimensions)).sum( + dim=kernel_dim) + stabilizer) + else: + return torch.exp(self.log_amp) * ((order_variances.unsqueeze(-1) * e_n.narrow(kernel_dim, 1, num_dimensions)).sum(dim=kernel_dim) + stabilizer) + + + # def grad(self, x1, x2=None): + # if x2 is None: + # x2 = x1 + # diffs = (x1[:, self.num_discrete:] - x2[:, self.num_discrete:])/self.lengthscales + # return diffs.t() + +if __name__ == '__main__': + pass diff --git a/GPmodel/kernels/mixedkernel.py b/GPmodel/kernels/mixedkernel.py new file mode 100644 index 0000000..6982e05 --- /dev/null +++ b/GPmodel/kernels/mixedkernel.py @@ -0,0 +1,33 @@ +import torch + +from GPmodel.modules.gp_modules import GPModule + + +class MixedKernel(GPModule): + + def __init__(self, log_order_variances, grouped_log_beta, fourier_freq_list, fourier_basis_list, lengthscales, num_discrete, num_continuous): + super(MixedKernel, self).__init__() + self.log_amp = torch.FloatTensor(1) + self.log_order_variances = log_order_variances # torch.ones(size=(num_discrete + num_continuous, )) # one for each combination of interaction + self.grouped_log_beta = grouped_log_beta + self.fourier_freq_list = fourier_freq_list + self.fourier_basis_list = fourier_basis_list + self.lengthscales = lengthscales + self.num_discrete = num_discrete + self.num_continuous = num_continuous + assert self.log_order_variances.size(0) == self.num_continuous + self.num_discrete, "order variances are not properly initialized" + assert self.lengthscales.size(0) == self.num_continuous, "lengthscales is not properly initialized" + assert self.grouped_log_beta.size(0) == self.num_discrete, "beta is not properly initialized" + + def n_params(self): + return 1 + + def param_to_vec(self): + return self.log_amp.clone() + + def vec_to_param(self, vec): + assert vec.numel() == 1 # self.num_discrete + self.num_continuous + self.log_amp = vec[:1].clone() + + def forward(self, input1, input2=None): + raise NotImplementedError diff --git a/GPmodel/likelihoods/__init__.py b/GPmodel/likelihoods/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/likelihoods/gaussian.py b/GPmodel/likelihoods/gaussian.py new file mode 100755 index 0000000..e4fd16b --- /dev/null +++ b/GPmodel/likelihoods/gaussian.py @@ -0,0 +1,31 @@ +import torch + +from GPmodel.likelihoods.likelihood import Likelihood + + +class GaussianLikelihood(Likelihood): + + def __init__(self): + super(GaussianLikelihood, self).__init__() + self.log_noise_var = torch.FloatTensor(1) + self.noise_scale = 0.1 + + def n_params(self): + return 1 + + def param_to_vec(self): + return self.log_noise_var.clone() + + def vec_to_param(self, vec): + self.log_noise_var = vec.clone() + + def forward(self, input): + return torch.exp(self.log_noise_var).repeat(input.size(0)) + + def __repr__(self): + return self.__class__.__name__ + + +if __name__ == '__main__': + likelihood = GaussianLikelihood() + print(list(likelihood.parameters())) diff --git a/GPmodel/likelihoods/likelihood.py b/GPmodel/likelihoods/likelihood.py new file mode 100755 index 0000000..5b37175 --- /dev/null +++ b/GPmodel/likelihoods/likelihood.py @@ -0,0 +1,7 @@ +from GPmodel.modules.gp_modules import GPModule + + +class Likelihood(GPModule): + + def __init__(self): + super(Likelihood, self).__init__() diff --git a/GPmodel/means/__init__.py b/GPmodel/means/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/means/constant.py b/GPmodel/means/constant.py new file mode 100755 index 0000000..d265d25 --- /dev/null +++ b/GPmodel/means/constant.py @@ -0,0 +1,31 @@ +import torch + +from GPmodel.means.mean import Mean + + +class ConstantMean(Mean): + + def __init__(self): + super(ConstantMean, self).__init__() + self.const_mean = torch.FloatTensor(1) + + def n_params(self): + return 1 + + def param_to_vec(self): + return self.const_mean.clone() + + def vec_to_param(self, vec): + self.const_mean = vec.clone() + + def forward(self, input): + # print("input ", input) + return self.const_mean * input.new_ones(input.size(0), 1).float() + + def __repr__(self): + return self.__class__.__name__ + + +if __name__ == '__main__': + likelihood = ConstantMean() + print(list(likelihood.parameters())) diff --git a/GPmodel/means/mean.py b/GPmodel/means/mean.py new file mode 100755 index 0000000..a7cd742 --- /dev/null +++ b/GPmodel/means/mean.py @@ -0,0 +1,7 @@ +from GPmodel.modules.gp_modules import GPModule + + +class Mean(GPModule): + + def __init__(self): + super(Mean, self).__init__() diff --git a/GPmodel/models/__init__.py b/GPmodel/models/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/models/gp.py b/GPmodel/models/gp.py new file mode 100755 index 0000000..c00154d --- /dev/null +++ b/GPmodel/models/gp.py @@ -0,0 +1,34 @@ +import torch +from GPmodel.modules.gp_modules import GPModule + + +class GP(GPModule): + def __init__(self, **kwargs): + super(GP, self).__init__() + + def init_param(self, output_data): + raise NotImplementedError + + def n_params(self): + cnt = 0 + print("heree") + for param in self.parameters(): + cnt += param.numel() + return cnt + + def param_to_vec(self): + flat_param_list = [] + for m in self.children(): + # print("************* [children] ************") + # print(m, m.param_to_vec()) + flat_param_list.append(m.param_to_vec()) + return torch.cat(flat_param_list) + + def vec_to_param(self, vec): + # print("vec", vec) + ind = 0 + for m in self.children(): + jump = m.n_params() + # print("jump: ", jump) + m.vec_to_param(vec[ind:ind+jump]) + ind += jump diff --git a/GPmodel/models/gp_regression.py b/GPmodel/models/gp_regression.py new file mode 100755 index 0000000..3a90d53 --- /dev/null +++ b/GPmodel/models/gp_regression.py @@ -0,0 +1,21 @@ +import torch + +from GPmodel.likelihoods.gaussian import GaussianLikelihood +from GPmodel.means.constant import ConstantMean +from GPmodel.models.gp import GP + + +class GPRegression(GP): + + def __init__(self, kernel, mean=ConstantMean()): + super(GPRegression, self).__init__() + self.kernel = kernel + self.mean = mean + self.likelihood = GaussianLikelihood() + + def init_param(self, output_data): + output_mean = torch.mean(output_data).item() + output_log_var = (0.5 * torch.var(output_data)).log().item() + self.kernel.log_amp.fill_(output_log_var) + self.mean.const_mean.fill_(output_mean) + self.likelihood.log_noise_var.fill_(output_mean / 1000.0) diff --git a/GPmodel/modules/__init__.py b/GPmodel/modules/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/modules/gp_modules.py b/GPmodel/modules/gp_modules.py new file mode 100755 index 0000000..4fc216d --- /dev/null +++ b/GPmodel/modules/gp_modules.py @@ -0,0 +1,19 @@ +from torch.nn.modules.module import Module + + +class GPModule(Module): + + def __init__(self): + super(GPModule, self).__init__() + + def n_params(self): + raise NotImplementedError + + def param_to_vec(self): + raise NotImplementedError + + def vec_to_param(self, vec): + raise NotImplementedError + + def __repr__(self): + return self.__class__.__name__ diff --git a/GPmodel/sampler/__init__.py b/GPmodel/sampler/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/GPmodel/sampler/priors.py b/GPmodel/sampler/priors.py new file mode 100755 index 0000000..196c75d --- /dev/null +++ b/GPmodel/sampler/priors.py @@ -0,0 +1,113 @@ +import numpy as np +from scipy.special import gammaln + +from GPmodel.sampler.tool_partition import compute_group_size + +# For numerical stability in exponential +LOG_LOWER_BND = -12.0 +LOG_UPPER_BND = 20.0 +# For sampling stability +STABLE_MEAN_RNG = 1.0 +# Hyperparameter for graph factorization +GRAPH_SIZE_LIMIT = 1024 + 2 + + +def log_prior_constmean(constmean, output_min, output_max): + """ + :param constmean: numeric(float) + :param output_min: numeric(float) + :param output_max: numeric(float) + :return: + """ + output_mid = (output_min + output_max) / 2.0 + output_rad = (output_max - output_min) * STABLE_MEAN_RNG / 2.0 + # Unstable parameter in sampling + if constmean < output_mid - output_rad or output_mid + output_rad < constmean: + return -float('inf') + # Uniform prior + # return 0 + # Truncated Gaussian + stable_dev = output_rad / 2.0 + return -np.log(stable_dev) - 0.5 * (constmean - output_mid) ** 2 / stable_dev ** 2 + + +def log_prior_noisevar(log_noise_var): + if log_noise_var < LOG_LOWER_BND or min(LOG_UPPER_BND, np.log(10000.0)) < log_noise_var: + return -float('inf') + return np.log(np.log(1.0 + (0.1 / np.exp(log_noise_var)) ** 2)) + + +def log_prior_kernelamp(log_amp, output_var, kernel_min, kernel_max): + """ + + :param log_amp: + :param output_var: numeric(float) + :param kernel_min: numeric(float) + :param kernel_max: numeric(float) + :return: + """ + if log_amp < LOG_LOWER_BND or min(LOG_UPPER_BND, np.log(10000.0)) < log_amp: + return -float('inf') + log_amp_lower = np.log(output_var) - np.log(kernel_max) + log_amp_upper = np.log(output_var) - np.log(max(kernel_min, 1e-100)) + log_amp_mid = 0.5 * (log_amp_upper + log_amp_lower) + log_amp_rad = 0.5 * (log_amp_upper - log_amp_lower) + log_amp_std = log_amp_rad / 2.0 + return -np.log(log_amp_std) - 0.5 * (log_amp - log_amp_mid) ** 2 / log_amp_std ** 2 + # return + # Uniform + # return 0 if kernel_min < output_var / amp < kernel_max else -float('inf') + # Gamma + # shape = output_var + # rate = 1.0 + # return shape * np.log(rate) - gammaln(shape) + (shape - 1.0) * log_amp - rate * np.exp(log_amp) + + +def log_prior_edgeweight(log_beta_i): + """ + + :param log_beta_i: numeric(float), ind-th element of 'log_beta' + :param dim: + :return: + """ + if log_beta_i < LOG_LOWER_BND or min(LOG_UPPER_BND, np.log(100.0)) < log_beta_i: + return -float('inf') + ## Gamma prior For sparsity-inducing, shape should be 1 + ## The higher the rate, the more sparsity is induced. + # shape = 1.0 + # rate = 3.0 + # return shape * np.log(rate) - gammaln(shape) + (shape - 1.0) * log_beta_i - rate * np.exp(log_beta_i) + + ## Horseshoe prior + tau = 5.0 + return np.log(np.log(1.0 + 2.0 / (np.exp(log_beta_i) / tau) ** 2)) + + ## Laplace prior + # scale = 0.5 + # return -np.exp(log_beta_i) / scale + + +def log_prior_partition(sorted_partition, n_vertices): + """ + Log of unnormalized density of given partition + this prior prefers well-spread partition, which is quantified by induced entropy. + Density is proportional to the entropy of a unnormalized probability vector consisting of [log(n_vertices in subgraph_i)]_i=1...N + :param sorted_partition: + :param n_vertices: 1D np.array + :return: + """ + if len(sorted_partition) == 1 or compute_group_size(sorted_partition=sorted_partition, n_vertices=n_vertices) > GRAPH_SIZE_LIMIT: + return -float('inf') + else: + prob_mass = np.array([np.sum(np.log(n_vertices[subset])) for subset in sorted_partition]) + prob_mass /= np.sum(prob_mass) + entropy_mass = -np.sum(prob_mass * np.log(prob_mass)) + max_log = np.sum(np.log(n_vertices)) + thr_log = np.log(GRAPH_SIZE_LIMIT) + n_chunk = int(np.floor(max_log / thr_log)) + prob_base = np.array([np.log(GRAPH_SIZE_LIMIT) for _ in range(n_chunk)] + [max_log - n_chunk * thr_log]) + prob_base /= np.sum(prob_base) + entropy_base = -np.sum(prob_base * np.log(prob_base)) + return np.log(entropy_mass - entropy_base) * 5 + # return np.log(entropy_mass) * 5 + diff --git a/GPmodel/sampler/sample_edgeweight.py b/GPmodel/sampler/sample_edgeweight.py new file mode 100755 index 0000000..d484e27 --- /dev/null +++ b/GPmodel/sampler/sample_edgeweight.py @@ -0,0 +1,62 @@ +import numpy as np + +import torch + +from GPmodel.inference.inference import Inference +from GPmodel.sampler.tool_partition import group_input +from GPmodel.sampler.tool_slice_sampling import univariate_slice_sampling +from GPmodel.sampler.priors import log_prior_edgeweight + + +def slice_edgeweight(model, input_data, output_data, n_vertices, log_beta, + sorted_partition, fourier_freq_list, fourier_basis_list, ind): + """ + Slice sampling the edgeweight(exp('log_beta')) at 'ind' in 'log_beta' vector + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param model: + :param input_data: + :param output_data: + :param n_vertices: 1d np.array + :param log_beta: + :param sorted_partition: Partition of {0, ..., K-1}, list of subsets(list) + :param fourier_freq_list: + :param fourier_basis_list: + :param ind: + :return: + """ + updated_subset_ind = [(ind in subset) for subset in sorted_partition].index(True) + updated_subset = sorted_partition[updated_subset_ind] + log_beta_rest = torch.sum(log_beta[updated_subset]) - log_beta[ind] + grouped_log_beta = torch.stack([torch.sum(log_beta[subset]) for subset in sorted_partition]) + model.kernel.grouped_log_beta = grouped_log_beta + model.kernel.fourier_freq_list = fourier_freq_list + model.kernel.fourier_basis_list = fourier_basis_list + grouped_input_data = group_input(input_data=input_data[:, :model.kernel.num_discrete], sorted_partition=sorted_partition, n_vertices=n_vertices)# need to understand this? + grouped_input_data = torch.cat((grouped_input_data, input_data[:, model.kernel.num_discrete:]), dim=1) + inference = Inference(train_data=(grouped_input_data, output_data), model=model) + + def logp(log_beta_i): + """ + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param log_beta_i: numeric(float) + :return: numeric(float) + """ + #if (log_beta_i < 0 or log_beta_i > 2): + # return -np.inf + log_prior = log_prior_edgeweight(log_beta_i) + if np.isinf(log_prior): + return log_prior + model.kernel.grouped_log_beta[updated_subset_ind] = log_beta_rest + log_beta_i + log_likelihood = float(-inference.negative_log_likelihood(hyper=model.param_to_vec())) + #return log_likelihood + return log_prior + log_likelihood + + x0 = float(log_beta[ind]) + x1 = univariate_slice_sampling(logp, x0) + log_beta[ind] = x1 + model.kernel.grouped_log_beta[updated_subset_ind] = log_beta_rest + x1 + return log_beta + + +if __name__ == '__main__': + pass diff --git a/GPmodel/sampler/sample_hyper.py b/GPmodel/sampler/sample_hyper.py new file mode 100755 index 0000000..a8dc423 --- /dev/null +++ b/GPmodel/sampler/sample_hyper.py @@ -0,0 +1,118 @@ +import numpy as np + +import torch + +from GPmodel.inference.inference import Inference +from GPmodel.sampler.tool_partition import group_input +from GPmodel.sampler.tool_slice_sampling import univariate_slice_sampling +from GPmodel.sampler.priors import log_prior_constmean, log_prior_noisevar, log_prior_kernelamp + + +def slice_hyper(model, input_data, output_data, n_vertices, sorted_partition): + """ + + :param model: + :param input_data: + :param output_data: + :return: + """ + grouped_input_data = group_input(input_data=input_data[:, :model.kernel.num_discrete], sorted_partition=sorted_partition, n_vertices=n_vertices)# need to understand this? + grouped_input_data = torch.cat((grouped_input_data, input_data[:, model.kernel.num_discrete:]), dim=1) + inference = Inference(train_data=(grouped_input_data, output_data), model=model) + #print("############# [slicing constmean] #############") + slice_constmean(inference) + #print("############# [slicing kernelamp] #############") + slice_kernelamp(inference) + #print("############# [slicing noise_var] #############") + slice_noisevar(inference) + + +def slice_constmean(inference): + """ + Slice sampling const_mean, this function does not need to return a sampled value + This directly modifies parameters in the argument 'inference.model.mean.const_mean' + :param inference: + :return: + """ + output_min = torch.min(inference.train_y).item() + output_max = torch.max(inference.train_y).item() + def logp(constmean): + """ + :param constmean: numeric(float) + :return: numeric(float) + """ + log_prior = log_prior_constmean(constmean, output_min=output_min, output_max=output_max) + if np.isinf(log_prior): + return log_prior + inference.model.mean.const_mean.fill_(constmean) + #print("data") + #print(inference.train_x) + log_likelihood = float(-inference.negative_log_likelihood(hyper=inference.model.param_to_vec())) + #print("&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& [here] &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&") + #print("!!!!!!!!!!!!!!! model to vec!!!!!!!!!!!!!!!!!") + # print(inference.model.param_to_vec()) + #print("log_prior:", log_prior) + #print("log_likelihood:", log_likelihood) + return log_prior + log_likelihood + + x0 = float(inference.model.mean.const_mean) + #print("&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& [here] &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&") + x1 = univariate_slice_sampling(logp, x0) + inference.model.mean.const_mean.fill_(x1) + return + + +def slice_noisevar(inference): + """ + Slice sampling log_noise_var, this function does not need to return a sampled value + This directly modifies parameters in the argument 'inference.model.likelihood.log_noise_var' + :param inference: + :return: + """ + + def logp(log_noise_var): + """ + :param log_noise_var: numeric(float) + :return: numeric(float) + """ + log_prior = log_prior_noisevar(log_noise_var) + if np.isinf(log_prior): + return log_prior + inference.model.likelihood.log_noise_var.fill_(log_noise_var) + log_likelihood = float(-inference.negative_log_likelihood(hyper=inference.model.param_to_vec())) + return log_prior + log_likelihood + + x0 = float(inference.model.likelihood.log_noise_var) + x1 = univariate_slice_sampling(logp, x0) + inference.model.likelihood.log_noise_var.fill_(x1) + return + + +def slice_kernelamp(inference): + """ + Slice sampling log_amp, this function does not need to return a sampled value + This directly modifies parameters in the argument 'inference.model.kernel.log_amp' + :param inference: + :return: + """ + output_var = torch.var(inference.train_y).item() + kernel_min = np.prod([torch.mean(torch.exp(-fourier_freq[-1])).item() / torch.mean(torch.exp(-fourier_freq)).item() for fourier_freq in inference.model.kernel.fourier_freq_list]) + kernel_max = np.prod([torch.mean(torch.exp(-fourier_freq[0])).item() / torch.mean(torch.exp(-fourier_freq)).item() for fourier_freq in inference.model.kernel.fourier_freq_list]) + #print(f"kernel_min : {kernel_min}, kernel_max: {kernel_max}") + def logp(log_amp): + """ + :param log_amp: numeric(float) + :return: numeric(float) + """ + log_prior = log_prior_kernelamp(log_amp, output_var, kernel_min, kernel_max) + #print(f"log_amp prior {log_prior}") + if np.isinf(log_prior): + return log_prior + inference.model.kernel.log_amp.fill_(log_amp) + log_likelihood = float(-inference.negative_log_likelihood(hyper=inference.model.param_to_vec())) + return log_prior + log_likelihood + + x0 = float(inference.model.kernel.log_amp) + x1 = univariate_slice_sampling(logp, x0) + inference.model.kernel.log_amp.fill_(x1) + return diff --git a/GPmodel/sampler/sample_mixed_posterior.py b/GPmodel/sampler/sample_mixed_posterior.py new file mode 100644 index 0000000..f167ace --- /dev/null +++ b/GPmodel/sampler/sample_mixed_posterior.py @@ -0,0 +1,102 @@ +import sys +import time +import copy +import numpy as np + +from GPmodel.sampler.sample_hyper import slice_hyper +from GPmodel.sampler.sample_edgeweight import slice_edgeweight +from GPmodel.sampler.tool_partition import direct_product +from config import PROGRESS_BAR_LEN +from GPmodel.sampler.slice_lengthscale import slice_lengthscale +from GPmodel.sampler.slice_log_order_variance import slice_log_order_variance + + +def posterior_sampling(model, input_data, output_data, n_vertices, adj_mat_list, log_order_variance, + log_beta, lengthscale, sorted_partition, n_sample, n_burn=0, n_thin=1): + """ + + :param model: + :param input_data: + :param output_data: + :param n_vertices: + :param adj_mat_list: + :param log_beta: + :param sorted_partition: + :param n_sample: + :param n_burn: + :param n_thin: + :return: + """ + + hyper_samples = [] + log_beta_samples = [] + log_order_variance_samples = [] + lengthscale_samples = [] + partition_samples = [] + freq_samples = [] + basis_samples = [] + edge_mat_samples = [] + + partition_sample = sorted_partition + log_beta_sample = log_beta + log_order_variance_sample = log_order_variance + lengthscale_sample = lengthscale + fourier_freq_list = model.kernel.fourier_freq_list + fourier_basis_list = model.kernel.fourier_basis_list + edge_mat_list = [direct_product(adj_mat_list, subset) for subset in sorted_partition] # kronecker operation + + n_sample_total = n_sample * n_thin + n_burn + n_digit = int(np.ceil(np.log(n_sample_total) / np.log(10))) + for s in range(0, n_sample_total): + #print("*************[ slicing hyper parameter ]*************") + slice_hyper(model, input_data, output_data, n_vertices, sorted_partition=partition_sample) + #print("model.mean:", model.mean.const_mean) + #print("model likelihood:", model.likelihood.log_noise_var) + #print("model kernel:", model.kernel.log_amp) + # In 'Batched High-dimensional Bayesian Optimization via Structural Kernel Learning', + # similar additive structure is updated for every 50 iterations(evaluations) + # This may be due to too much variability if decomposition is learned every iterations. + # Thus, when multiple points are sampled, sweeping all inds for each sample may be not a good practice + # For example if there are 50 variables and 10 samples are needed, + # then after shuffling indices, and 50/10 thinning can be used. + + shuffled_beta_ind = list(range(len(adj_mat_list))) + np.random.shuffle(shuffled_beta_ind) + for beta_ind in shuffled_beta_ind: + # In each sampler, model.kernel fourier_freq_list, fourier_basis_list are updated. + log_beta_sample = slice_edgeweight(model, input_data, output_data, n_vertices, + log_beta=log_beta_sample, sorted_partition=partition_sample, + fourier_freq_list=fourier_freq_list, fourier_basis_list=fourier_basis_list, + ind=beta_ind) + + shuffled_ls_ind = list(range(model.kernel.num_continuous)) + np.random.shuffle(shuffled_ls_ind) + for ls_ind in shuffled_ls_ind: + lengthscale_sample = slice_lengthscale(model, input_data, output_data, n_vertices, + lengthscale_sample=lengthscale_sample, sorted_partition=partition_sample, + fourier_freq_list=fourier_freq_list, fourier_basis_list=fourier_basis_list, + ind=ls_ind) + shuffled_lov_ind = list(range(model.kernel.num_continuous + model.kernel.num_discrete)) + np.random.shuffle(shuffled_lov_ind) + for lov_ind in shuffled_lov_ind: + log_order_variance_sample = slice_log_order_variance(model, input_data, output_data, n_vertices, + log_order_variance=log_order_variance_sample, sorted_partition=partition_sample, + fourier_freq_list=fourier_freq_list, fourier_basis_list=fourier_basis_list, + ind=lov_ind) + if s >= n_burn and (s - n_burn + 1) % n_thin == 0: + hyper_samples.append(model.param_to_vec()) + log_beta_samples.append(log_beta_sample.clone()) + lengthscale_samples.append(lengthscale_sample.clone()) + log_order_variance_samples.append(log_order_variance_sample.clone()) + partition_samples.append(copy.deepcopy(partition_sample)) + freq_samples.append([elm.clone() for elm in fourier_freq_list]) + basis_samples.append([elm.clone() for elm in fourier_basis_list]) + edge_mat_samples.append([elm.clone() for elm in edge_mat_list]) + progress_mark_len = int((s + 1.0) / n_sample_total * PROGRESS_BAR_LEN) + fmt_str = '(%s) %3d%% (%' + str(n_digit) + 'd of %d) |' \ + + '#' * progress_mark_len + '-' * (PROGRESS_BAR_LEN - progress_mark_len) + '|' + progress_str = fmt_str % (time.strftime('%H:%M:%S', time.gmtime()), + int((s + 1.0) / n_sample_total * 100), s + 1, n_sample_total) + sys.stdout.write(('\b' * len(progress_str)) + progress_str) + sys.stdout.flush() + return hyper_samples, log_order_variance_samples, log_beta_samples, lengthscale_samples, partition_samples, freq_samples, basis_samples, edge_mat_samples diff --git a/GPmodel/sampler/sample_posterior.py b/GPmodel/sampler/sample_posterior.py new file mode 100755 index 0000000..062cadb --- /dev/null +++ b/GPmodel/sampler/sample_posterior.py @@ -0,0 +1,88 @@ +import sys +import time +import copy +import numpy as np + +from graphGP.sampler.sample_hyper import slice_hyper +from graphGP.sampler.sample_edgeweight import slice_edgeweight +from graphGP.sampler.tool_partition import direct_product +from config import PROGRESS_BAR_LEN + + +def posterior_sampling(model, input_data, output_data, n_vertices, adj_mat_list, + log_beta, sorted_partition, n_sample, n_burn=0, n_thin=1): + """ + + :param model: + :param input_data: + :param output_data: + :param n_vertices: + :param adj_mat_list: + :param log_beta: + :param sorted_partition: + :param n_sample: + :param n_burn: + :param n_thin: + :return: + """ + hyper_samples = [] + log_beta_samples = [] + log_lengthscale_samples = [] + partition_samples = [] + freq_samples = [] + basis_samples = [] + edge_mat_samples = [] + + partition_sample = sorted_partition + log_beta_sample = log_beta + fourier_freq_list = model.kernel.fourier_freq_list + fourier_basis_list = model.kernel.fourier_basis_list + edge_mat_list = [direct_product(adj_mat_list, subset) for subset in sorted_partition] # kronecker operation + + n_sample_total = n_sample * n_thin + n_burn + n_digit = int(np.ceil(np.log(n_sample_total) / np.log(10))) + for s in range(0, n_sample_total): + slice_hyper(model, input_data, output_data, n_vertices, sorted_partition=partition_sample) + # In 'Batched High-dimensional Bayesian Optimization via Structural Kernel Learning', + # similar additive structure is updated for every 50 iterations(evaluations) + # This may be due to too much variability if decomposition is learned every iterations. + # Thus, when multiple points are sampled, sweeping all inds for each sample may be not a good practice + # For example if there are 50 variables and 10 samples are needed, + # then after shuffling indices, and 50/10 thinning can be used. + + shuffled_beta_ind = list(range(len(adj_mat_list))) + np.random.shuffle(shuffled_beta_ind) + for beta_ind in shuffled_beta_ind: + # In each sampler, model.kernel fourier_freq_list, fourier_basis_list are updated. + log_beta_sample = slice_edgeweight(model, input_data, output_data, n_vertices, + log_beta=log_beta_sample, sorted_partition=partition_sample, + fourier_freq_list=fourier_freq_list, fourier_basis_list=fourier_basis_list, + ind=beta_ind) + + shuffled_ls_ind = list(range(model.kernel.num_continuous)) + np.random.shuffle(shuffled_ls_ind) + for ls_ind in shuffled_ls_ind: + # In each sampler, model.kernel fourier_freq_list, fourier_basis_list are updated. + log_lengthscale_sample = slice_edgeweight(model, input_data, output_data, n_vertices, + log_lengthscale_sample=log_lengthscale_sample, sorted_partition=partition_sample, + fourier_freq_list=fourier_freq_list, fourier_basis_list=fourier_basis_list, + ind=ls_ind) + + if s >= n_burn and (s - n_burn + 1) % n_thin == 0: + hyper_samples.append(model.param_to_vec()) + print("model param to vec") + print(model.param_to_vec()) + log_beta_samples.append(log_beta_sample.clone()) + log_lengthscale_sample.append(log_ls_sample.clone()) + partition_samples.append(copy.deepcopy(partition_sample)) + freq_samples.append([elm.clone() for elm in fourier_freq_list]) + basis_samples.append([elm.clone() for elm in fourier_basis_list]) + edge_mat_samples.append([elm.clone() for elm in edge_mat_list]) + progress_mark_len = int((s + 1.0) / n_sample_total * PROGRESS_BAR_LEN) + fmt_str = '(%s) %3d%% (%' + str(n_digit) + 'd of %d) |' \ + + '#' * progress_mark_len + '-' * (PROGRESS_BAR_LEN - progress_mark_len) + '|' + progress_str = fmt_str % (time.strftime('%H:%M:%S', time.gmtime()), + int((s + 1.0) / n_sample_total * 100), s + 1, n_sample_total) + sys.stdout.write(('\b' * len(progress_str)) + progress_str) + sys.stdout.flush() + return hyper_samples, log_beta_samples, partition_samples, freq_samples, basis_samples, edge_mat_samples diff --git a/GPmodel/sampler/slice_lengthscale.py b/GPmodel/sampler/slice_lengthscale.py new file mode 100644 index 0000000..4f73ef8 --- /dev/null +++ b/GPmodel/sampler/slice_lengthscale.py @@ -0,0 +1,54 @@ +import numpy as np + +import torch + +from GPmodel.inference.inference import Inference +from GPmodel.sampler.tool_partition import group_input +from GPmodel.sampler.tool_slice_sampling import univariate_slice_sampling +from GPmodel.sampler.priors import log_prior_edgeweight + +# no prior as such but best to keep the values within 0 and a maximum (100) + + +def slice_lengthscale(model, input_data, output_data, n_vertices, lengthscale_sample, + sorted_partition, fourier_freq_list, fourier_basis_list, ind): + """ + Slice sampling at 'ind' in 'log_lengthscale' vector + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param model: + :param input_data: + :param output_data: + :param n_vertices: 1d np.array + :param log_beta: + :param sorted_partition: Partition of {0, ..., K-1}, list of subsets(list) + :param fourier_freq_list: + :param fourier_basis_list: + :param ind: + :return: + """ + grouped_input_data = group_input(input_data=input_data[:, :model.kernel.num_discrete], sorted_partition=sorted_partition, n_vertices=n_vertices)# need to understand this? + grouped_input_data = torch.cat((grouped_input_data, input_data[:, model.kernel.num_discrete:]), dim=1) + inference = Inference(train_data=(grouped_input_data, output_data), model=model) + + def logp(ls): + """ + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param log_beta_i: numeric(float) + :return: numeric(float) + """ + if (ls < 0 or ls > 2): + return -np.inf#float('-inf') + inference.model.kernel.lengthscales[ind] = ls + log_likelihood = float(-inference.negative_log_likelihood(hyper=model.param_to_vec())) + return log_likelihood + + x0 = float(lengthscale_sample[ind]) + x1 = univariate_slice_sampling(logp, x0) + lengthscale_sample[ind] = x1 + model.kernel.lengthscales[ind] = x1 + return lengthscale_sample + + + +if __name__ == '__main__': + pass diff --git a/GPmodel/sampler/slice_log_order_variance.py b/GPmodel/sampler/slice_log_order_variance.py new file mode 100644 index 0000000..b05bf47 --- /dev/null +++ b/GPmodel/sampler/slice_log_order_variance.py @@ -0,0 +1,52 @@ +import numpy as np + +import torch + +from GPmodel.inference.inference import Inference +from GPmodel.sampler.tool_partition import group_input +from GPmodel.sampler.tool_slice_sampling import univariate_slice_sampling +from GPmodel.sampler.priors import log_prior_edgeweight + + +def slice_log_order_variance(model, input_data, output_data, n_vertices, log_order_variance, + sorted_partition, fourier_freq_list, fourier_basis_list, ind): + """ + Slice sampling the edgeweight(exp('log_beta')) at 'ind' in 'log_beta' vector + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param model: + :param input_data: + :param output_data: + :param n_vertices: 1d np.array + :param log_beta: + :param sorted_partition: Partition of {0, ..., K-1}, list of subsets(list) + :param fourier_freq_list: + :param fourier_basis_list: + :param ind: + :return: + """ + grouped_input_data = group_input(input_data=input_data[:, :model.kernel.num_discrete], sorted_partition=sorted_partition, n_vertices=n_vertices)# need to understand this? + grouped_input_data = torch.cat((grouped_input_data, input_data[:, model.kernel.num_discrete:]), dim=1) + inference = Inference(train_data=(grouped_input_data, output_data), model=model) + + def logp(log_order_variance_i): + """ + Note that model.kernel members (fourier_freq_list, fourier_basis_list) are updated. + :param log_beta_i: numeric(float) + :return: numeric(float) + """ + log_prior = log_prior_edgeweight(log_order_variance_i) + if np.isinf(log_prior): + return log_prior + model.kernel.log_order_variances[ind] = log_order_variance_i + log_likelihood = float(-inference.negative_log_likelihood(hyper=model.param_to_vec())) + return log_prior + log_likelihood + + x0 = float(log_order_variance[ind]) + x1 = univariate_slice_sampling(logp, x0) + log_order_variance[ind] = x1 + model.kernel.log_order_variances[ind] = x1 + return log_order_variance + + +if __name__ == '__main__': + pass diff --git a/GPmodel/sampler/tool_partition.py b/GPmodel/sampler/tool_partition.py new file mode 100755 index 0000000..5a775b2 --- /dev/null +++ b/GPmodel/sampler/tool_partition.py @@ -0,0 +1,123 @@ +import numpy as np + +import torch + + +def np_kron(mat1, mat2): + """ + In order to check the function below 'def kronecker', numpy kron is used + :param mat1: 2d torch.Tensor + :param mat2: 2d torch.Tensor + :return: kronecker product of mat1 and mat2 + """ + np_mat1 = mat1.numpy() + np_mat2 = mat2.numpy() + return torch.from_numpy(np.kron(np_mat1, np_mat2)) + + +def kronecker(mat1, mat2): + """ + kronecker product between 2 2D tensors + :param mat1: 2d torch.Tensor + :param mat2: 2d torch.Tensor + :return: kronecker product of mat1 and mat2 + """ + s1 = mat1.size() + s2 = mat2.size() + return torch.ger(mat1.view(-1), mat2.view(-1)).reshape(*(s1 + s2)).permute([0, 2, 1, 3]).reshape(s1[0] * s2[0], s1[1] * s2[1]) + + +def sort_partition(partition): + """ + given partition is ordered to have unique representation + :param partition: list of subsets + :return: each subset is represented as an ordered list, subsets are ordered according to their smallest elements + """ + lowest_ind_key_dict = {min(subset): sorted(subset) for subset in partition} + sorted_partition = [] + for k in sorted(lowest_ind_key_dict.keys()): + sorted_partition.append(lowest_ind_key_dict[k]) + return sorted_partition + + +def compute_unit_in_group(sorted_partition, n_vertices): + """ + In order to convert between grouped variable and original ungrouped variable, units are necessary + e.g) when C1, C2, C5 are grouped and C1, C2, C5 have 3, 4, 5 n_vertices respectively, + then (c1, c2, c5) in ungrouped value is equivalent to c1 x 4 x 5 + c2 x 5 + c3 in grouped value, + here 4 x 5 is the unit for c1, 5 is unit for c2, 1 is unit for c3 + :param sorted_partition: list of subsets, elements in a subset are ordered, subsets are ordered by their smallest elements + :param n_vertices: 1d np.array + :return: unit is given as the same order corresponding to sorted_partition + """ + unit_in_group = [] + for subset in sorted_partition: + ind_units = list(np.flip(np.cumprod((n_vertices[subset][1:][::-1])))) + [1] + unit_in_group.append(ind_units) + return unit_in_group + + +def compute_group_size(sorted_partition, n_vertices): + """ + Return the size of the largest subgraph (a product of strong product) + :param sorted_partition: + :param n_vertices: 1d np.array + :return: + """ + complexity = sum([np.prod(n_vertices[subset]) ** 3 for subset in sorted_partition]) + max_size = max([np.prod(n_vertices[subset]) for subset in sorted_partition]) + return max_size + + +def group_input(input_data, sorted_partition, n_vertices): + """ + + :param input_data: 2D torch.Tensor, size(0) : number of data, size(1) : number of original variables + :param sorted_partition: list of subsets, elements in each subset are ordered, subsets are ordered by their smallest elements + :param unit_in_group: compute_unit_in_group(sorted_partition, n_vertices) + :return: 2D torch.Tensor, size(0) : number of data, size(1) : number of grouped variables + """ + unit_in_group = compute_unit_in_group(sorted_partition, n_vertices) + grouped_input = input_data.new_zeros(input_data.size(0), len(sorted_partition)) + for g in range(len(sorted_partition)): + for ind, unit in zip(sorted_partition[g], unit_in_group[g]): + grouped_input[:, g] += input_data[:, ind] * unit + return grouped_input + + +def ungroup_input(grouped_input, sorted_partition, n_vertices): + """ + + :param grouped_input: 2D torch.Tensor, size(0) : number of data, size(1) : number of grouped variables + :param sorted_partition: list of subsets, elements in each subset are ordered, subsets are ordered by their smallest elements + :param unit_in_group: compute_unit_in_group(sorted_partition, n_vertices) + :return: 2D torch.Tensor, size(0) : number of data, size(1) : number of original variables + """ + unit_in_group = compute_unit_in_group(sorted_partition, n_vertices) + input_data = grouped_input.new_zeros(grouped_input.size(0), sum([len(subset) for subset in sorted_partition])) + for g in range(len(sorted_partition)): + subset = sorted_partition[g] + unit = unit_in_group[g] + elm_ind = 0 + remainder = grouped_input[:, g] + while elm_ind < len(subset) - 1: + input_data[:, subset[elm_ind]] = remainder // unit[elm_ind] + remainder = remainder % unit[elm_ind] + elm_ind += 1 + input_data[:, subset[-1]] = remainder + return input_data + + +def direct_product(adj_mat_list, subset): + """ + Adjacency matrix of direct product + :param adj_mat_list: list of 2D tensors or adjacency matrices + :param subset: list of subsets, elements in each subset are ordered, subsets are ordered by their smallest elements + :return: + """ + elm = subset[0] + grouped_adjacency = adj_mat_list[elm] + for ind in range(1, len(subset)): + elm = subset[ind] + grouped_adjacency = kronecker(grouped_adjacency, adj_mat_list[elm]) + return grouped_adjacency diff --git a/GPmodel/sampler/tool_slice_sampling.py b/GPmodel/sampler/tool_slice_sampling.py new file mode 100755 index 0000000..e8a0a3b --- /dev/null +++ b/GPmodel/sampler/tool_slice_sampling.py @@ -0,0 +1,93 @@ +import numpy as np + + +def univariate_slice_sampling(logp, x0, width=1.0, max_steps_out=10): + """ + Univariate Slice Sampling using doubling scheme + :param logp: numeric(float) -> numeric(float), a log density function + :param x0: numeric(float) + :param width: + :param max_steps_out: + :return: numeric(float), sampled x1 + """ + # print("x0", x0) + for scaled_width in np.array([0.9, 0.8, 0.7, 0.6, 0.4, 0.3, 0.2, 0.1]) * width: + + lower = x0 - scaled_width * np.random.rand() + upper = lower + scaled_width + llh0 = logp(x0) + #print("llh0:", llh0) + slice_h = np.log(np.random.rand()) + llh0 + llh_record = {} + + # Step Out (doubling) + #print("Lower and upper") + #print(lower) + steps_out = 0 + logp_lower = logp(lower) + logp_upper = logp(upper) + llh_record[float(lower)] = logp_lower + llh_record[float(upper)] = logp_upper + while (logp_lower > slice_h or logp_upper > slice_h) and (steps_out < max_steps_out): + if np.random.rand() < 0.5: + lower -= (upper - lower) + else: + upper += (upper - lower) + steps_out += 1 + try: + logp_lower = llh_record[float(lower)] + except KeyError: + logp_lower = logp(lower) + llh_record[float(lower)] = logp_lower + try: + logp_upper = llh_record[float(upper)] + except KeyError: + logp_upper = logp(lower) + llh_record[float(upper)] = logp_upper + + # Shrinkage + start_upper = upper + start_lower = lower + n_steps_in = 0 + while not np.isclose(lower, upper): + x1 = (upper - lower) * np.random.rand() + lower + llh1 = logp(x1) + if llh1 > slice_h and accept(logp, x0, x1, slice_h, scaled_width, start_lower, start_upper, llh_record): + return x1 + else: + if x1 < x0: + lower = x1 + else: + upper = x1 + n_steps_in += 1 + # raise RuntimeError('Shrinkage collapsed to a degenerated interval(point)') + return x0 # just returning original value + + +def accept(logp, x0, x1, slice_h, width, lower, upper, llh_record): + acceptance = False + while upper - lower > 1.1 * width: + mid = (lower + upper) / 2.0 + if (x0 < mid and x1 >= mid) or (x0 >= mid and x1 < mid): + acceptance = True + if x1 < mid: + upper = mid + else: + lower = mid + try: + logp_lower = llh_record[float(lower)] + except KeyError: + logp_lower = logp(lower) + llh_record[float(lower)] = logp_lower + try: + logp_upper = llh_record[float(upper)] + except KeyError: + logp_upper = logp(upper) + llh_record[float(upper)] = logp_upper + if acceptance and slice_h >= logp_lower and slice_h >= logp_upper: + return False + return True + + +if __name__ == '__main__': + pass \ No newline at end of file diff --git a/acquisition/__init__.py b/acquisition/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/acquisition/acquisition_functions.py b/acquisition/acquisition_functions.py new file mode 100755 index 0000000..a1e6d36 --- /dev/null +++ b/acquisition/acquisition_functions.py @@ -0,0 +1,18 @@ +import torch +from torch.distributions.normal import Normal + + +def expected_improvement(mean, var, reference): + """ + expected_improvement for minimization problems + :param mean: + :param var: + :param reference: + :return: + """ + predictive_normal = Normal(mean.new_zeros(mean.size()), mean.new_ones(mean.size())) + std = torch.sqrt(var) + standardized = (-mean + reference) / std + return (std * torch.exp(predictive_normal.log_prob(standardized)) + (-mean + reference) * predictive_normal.cdf(standardized)).clamp(min=0) + + diff --git a/acquisition/acquisition_marginalization.py b/acquisition/acquisition_marginalization.py new file mode 100755 index 0000000..b88d486 --- /dev/null +++ b/acquisition/acquisition_marginalization.py @@ -0,0 +1,77 @@ +import torch + +from GPmodel.kernels.mixeddiffusionkernel import MixedDiffusionKernel +from GPmodel.models.gp_regression import GPRegression +from GPmodel.inference.inference import Inference +from GPmodel.sampler.tool_partition import group_input + +from acquisition.acquisition_functions import expected_improvement + + +def acquisition_expectation(x, inference_samples, partition_samples, n_vertices, acquisition_func=expected_improvement, + reference=None): + """ + Using posterior samples, the acquisition function is also averaged over posterior samples + :param x: 1d or 2d tensor + :param inference_samples: inference method for each posterior sample + :param partition_samples: + :param n_vertices: + :param acquisition_func: + :param reference: + :return: + """ + if x.dim() == 1: + x = x.unsqueeze(0) + + acquisition_sample_list = [] + for s in range(len(inference_samples)): + hyper = inference_samples[s].model.param_to_vec() + grouped_x = group_input(x[:, :inference_samples[s].model.kernel.num_discrete], sorted_partition=partition_samples[s], n_vertices=n_vertices) + grouped_x = torch.cat((grouped_x, x[:, inference_samples[s].model.kernel.num_discrete:]), dim=1) + pred_dist = inference_samples[s].predict(grouped_x, hyper=hyper, verbose=False) + pred_mean_sample = pred_dist[0].detach() + pred_var_sample = pred_dist[1].detach() + acquisition_sample_list.append(acquisition_func(pred_mean_sample[:, 0], pred_var_sample[:, 0], + reference=reference)) + + return torch.stack(acquisition_sample_list, 1).sum(1, keepdim=True) + + +def inference_sampling(input_data, output_data, n_vertices, hyper_samples, log_order_variance_samples, log_beta_samples, lengthscales_samples, partition_samples, + freq_samples, basis_samples, num_discrete, num_continuous): + """ + """ + inference_samples = [] + for s in range(len(hyper_samples)): + grouped_log_beta = torch.stack([torch.sum(log_beta_samples[s][subset]) for subset in partition_samples[s]]) + kernel = MixedDiffusionKernel(log_order_variances = log_order_variance_samples[s], grouped_log_beta=grouped_log_beta, fourier_freq_list=freq_samples[s], + fourier_basis_list=basis_samples[s], lengthscales=lengthscales_samples[s], + num_discrete=num_discrete, num_continuous=num_continuous) + model = GPRegression(kernel=kernel) + model.vec_to_param(hyper_samples[s]) + grouped_input_data = group_input(input_data=input_data[:, :model.kernel.num_discrete], sorted_partition=partition_samples[s], n_vertices=n_vertices)# need to understand this? + grouped_input_data = torch.cat((grouped_input_data, input_data[:, model.kernel.num_discrete:]), dim=1) + inference = Inference((grouped_input_data, output_data), model=model) + inference_samples.append(inference) + return inference_samples + + +def prediction_statistic(x, inference_samples, partition_samples, n_vertices): + if x.dim() == 1: + x = x.unsqueeze(0) + mean_sample_list = [] + std_sample_list = [] + var_sample_list = [] + for s in range(len(inference_samples)): + grouped_x = group_input(input_data=x[:, :inference_samples[s].model.kernel.num_discrete], sorted_partition=partition_samples[s], n_vertices=n_vertices) + grouped_x = torch.cat((grouped_x, x[:, inference_samples[s].model.kernel.num_discrete:]), dim=1) + pred_dist = inference_samples[s].predict(grouped_x) + pred_mean_sample = pred_dist[0] + pred_var_sample = pred_dist[1] + pred_std_sample = pred_var_sample ** 0.5 + mean_sample_list.append(pred_mean_sample.data) + std_sample_list.append(pred_std_sample.data) + var_sample_list.append(pred_var_sample.data) + return torch.cat(mean_sample_list, 1).mean(1, keepdim=True),\ + torch.cat(std_sample_list, 1).mean(1, keepdim=True),\ + torch.cat(var_sample_list, 1).mean(1, keepdim=True) diff --git a/acquisition/acquisition_optimization.py b/acquisition/acquisition_optimization.py new file mode 100755 index 0000000..2631ac0 --- /dev/null +++ b/acquisition/acquisition_optimization.py @@ -0,0 +1,114 @@ +import os +import sys +import time +import psutil + +import numpy as np + +import torch +import torch.multiprocessing as mp + +from acquisition.acquisition_optimizers.starting_points import optim_inits +from acquisition.acquisition_optimizers.greedy_ascent import greedy_ascent +from acquisition.acquisition_optimizers.continuous_optimizer import cma_es_optimizer +from acquisition.acquisition_functions import expected_improvement +from acquisition.acquisition_marginalization import prediction_statistic +from acquisition.acquisition_optimizers.graph_utils import neighbors + + +MAX_N_ASCENT = float('inf') +N_CPU = os.cpu_count() +N_AVAILABLE_CORE = min(10, N_CPU) +N_SA_RUN = 10 + + +def next_evaluation(objective, x_opt, input_data, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func=expected_improvement, reference=None, parallel=None): + id_digit = np.ceil(np.log(np.prod(n_vertices)) / np.log(10)) + id_unit = torch.from_numpy(np.cumprod(np.concatenate([np.ones(1), n_vertices[:-1]])).astype(np.int)) + fmt_str = '\t %5.2f (id:%' + str(id_digit) + 'd) ==> %5.2f (id:%' + str(id_digit) + 'd)' + + start_time = time.time() + print('(%s) Acquisition function optimization initial points selection began' + % (time.strftime('%H:%M:%S', time.localtime(start_time)))) + + x_inits, acq_inits = optim_inits(objective, x_opt, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func, reference) + n_inits = x_inits.size(0) + assert n_inits % 2 == 0 + + end_time = time.time() + elapsed_time = end_time - start_time + print('(%s) Acquisition function optimization initial points selection ended - %s' + % (time.strftime('%H:%M:%S', time.localtime(end_time)), time.strftime('%H:%M:%S', time.localtime(elapsed_time)))) + + start_time = time.time() + print('(%s) Acquisition function optimization with %2d inits' + % (time.strftime('%H:%M:%S', time.localtime(start_time)), x_inits.size(0))) + + ga_args_list = [(x_inits[i], inference_samples, partition_samples, edge_mat_samples, + n_vertices, acquisition_func, MAX_N_ASCENT, reference) for i in range(n_inits)] + ga_start_time = time.time() + sys.stdout.write(' Greedy Ascent began at %s ' % time.strftime('%H:%M:%S', time.localtime(ga_start_time))) + if parallel: + with mp.Pool(processes=min(n_inits, N_CPU // 3)) as pool: + ga_result = [] + process_started = [False] * n_inits + process_running = [False] * n_inits + process_index = 0 + while process_started.count(False) > 0: + cpu_usage = psutil.cpu_percent(0.25) + run_more = (100.0 - cpu_usage) * float(psutil.cpu_count()) > 100.0 * N_AVAILABLE_CORE + if run_more: + ga_result.append(pool.apply_async(greedy_ascent, args=ga_args_list[process_index])) + process_started[process_index] = True + process_running[process_index] = True + process_index += 1 + while [not res.ready() for res in ga_result].count(True) > 0: + time.sleep(1) + ga_return_values = [res.get() for res in ga_result] + else: + ga_return_values = [greedy_ascent(*(ga_args_list[i])) for i in range(n_inits)] + + ga_args_list = [(objective, ga_return_values[i][0], ga_return_values[i][1], inference_samples, partition_samples, + n_vertices, expected_improvement, reference) for i in range(len(ga_return_values))] + + ga_return_values = [cma_es_optimizer(*(ga_args_list[i])) for i in range(len(ga_args_list))] + ga_opt_vrt, ga_opt_acq = zip(*ga_return_values) + sys.stdout.write('and took %s\n' % time.strftime('%H:%M:%S', time.localtime(time.time() - ga_start_time))) + + opt_vrt = list(ga_opt_vrt[:]) + opt_acq = list(ga_opt_acq[:]) + + end_time = time.time() + elapsed_time = end_time - start_time + print('(%s) Acquisition function optimization ended %s' + % (time.strftime('%H:%M:%S', time.localtime(end_time)), time.strftime('%H:%M:%S', time.localtime(elapsed_time)))) + + # argsort sorts in ascending order so it is negated to have descending order + acq_sort_inds = np.argsort(-np.array(opt_acq)) + suggestion = None + for i in range(len(opt_vrt)): + ind = acq_sort_inds[i] + if not torch.all(opt_vrt[ind] == input_data, dim=1).any(): + suggestion = opt_vrt[ind] + break + if suggestion is None: + print("random suggestion!!!") + for i in range(len(opt_vrt)): + ind = acq_sort_inds[i] + nbds = neighbors(opt_vrt[ind][:inference_samples[0].model.kernel.num_discrete], partition_samples, edge_mat_samples, n_vertices, uniquely=True) + nbds = torch.cat((nbds, opt_vrt[ind][inference_samples[0].model.kernel.num_discrete:].unsqueeze(0).repeat(nbds.size(0), 1)), dim=1) + #nbds = neighbors(opt_vrt[ind], partition_samples, edge_mat_samples, n_vertices, uniquely=True) + for j in range(nbds.size(0)): + if not torch.all(nbds[j] == input_data, dim=1).any(): + suggestion = nbds[j] + break + if suggestion is not None: + break + if suggestion is None: + suggestion = torch.cat(tuple([torch.randint(low=0, high=int(n_v), size=(1, 1)) for n_v in n_vertices]), dim=1).long() + if (torch.all(suggestion == input_data, dim = 1).any()): + suggestion = objective.generate_random_points(1) + mean, std, var = prediction_statistic(suggestion, inference_samples, partition_samples, n_vertices) + return suggestion, mean, std, var diff --git a/acquisition/acquisition_optimizers/__init__.py b/acquisition/acquisition_optimizers/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/acquisition/acquisition_optimizers/continuous_optimizer.py b/acquisition/acquisition_optimizers/continuous_optimizer.py new file mode 100644 index 0000000..1815d32 --- /dev/null +++ b/acquisition/acquisition_optimizers/continuous_optimizer.py @@ -0,0 +1,51 @@ +import torch + +from acquisition.acquisition_functions import expected_improvement +from acquisition.acquisition_marginalization import acquisition_expectation +import numpy as np +import cma +import time +import scipy.optimize as spo +from functools import partial + + +def continuous_acquisition_expectation(x_continuous, discrete_part, inference_samples, partition_samples, + n_vertices, acquisition_func, reference, batch=False): + if batch: + eval_x = torch.from_numpy(np.concatenate((np.tile(discrete_part, (len(x_continuous), 1)), x_continuous), axis = 1)).float() + results = acquisition_expectation(eval_x,inference_samples, partition_samples, n_vertices,acquisition_func, reference) + return np.array(results) + else: + eval_x = torch.from_numpy(np.concatenate((discrete_part, x_continuous))).float() + print(acquisition_expectation(eval_x, inference_samples, partition_samples, n_vertices, + acquisition_func, reference)[0].numpy()) + return acquisition_expectation(eval_x, inference_samples, partition_samples, n_vertices, + acquisition_func, reference)[0].numpy() + + + +def cma_es_optimizer(objective, x_init, max_acquisition, inference_samples, partition_samples, n_vertices, acquisition_func=expected_improvement, reference=None): + cont_bounds = [objective.problem.lower_bounds[objective.num_discrete:], objective.problem.upper_bounds[objective.num_discrete:]] + start_time = time.time() + es = cma.CMAEvolutionStrategy(x0=x_init[objective.num_discrete:],sigma0=0.1,inopts={'bounds': cont_bounds, "popsize": 50},) + iter = 1 + total_time_in_acq = 0 + while not es.stop(): + iter += 1 + xs = es.ask() + X = torch.tensor(xs).float() + # evaluate the acquisition function (optimizer assumes we're minimizing) + temp_time = time.time() + Y = -1 * continuous_acquisition_expectation(xs, x_init[:objective.num_discrete].numpy(), + inference_samples, partition_samples, n_vertices, acquisition_func, reference, batch=True) + total_time_in_acq += time.time() - temp_time + es.tell(xs, Y) # return the result to the optimizer + if (iter > 10): + break + best_x = torch.from_numpy(es.best.x).float() + if -1*es.best.f > max_acquisition: + return torch.cat((x_init[:objective.num_discrete], best_x), dim=0), -1*es.best.f + else: + return x_init, max_acquisition + + diff --git a/acquisition/acquisition_optimizers/graph_utils.py b/acquisition/acquisition_optimizers/graph_utils.py new file mode 100755 index 0000000..6ddb2f3 --- /dev/null +++ b/acquisition/acquisition_optimizers/graph_utils.py @@ -0,0 +1,50 @@ +import torch + +from GPmodel.sampler.tool_partition import group_input, ungroup_input + + +def neighbors(x, partition_samples, edge_mat_samples, n_vertices, uniquely=False): + """ + + :param x: 1D Tensor + :param partition_samples: + :param edge_mat_samples: + :param n_vertices: + :param uniquely: + :return: + """ + nbds = x.new_empty((0, x.numel())) + for i in range(len(partition_samples)): + grouped_x = group_input(x.unsqueeze(0), partition_samples[i], n_vertices).squeeze(0) + grouped_nbd = _cartesian_neighbors(grouped_x, edge_mat_samples[i]) + nbd = ungroup_input(grouped_nbd, partition_samples[i], n_vertices) + added_ind = [] + if uniquely: + for j in range(nbd.size(0)): + if not torch.any(torch.all(nbds == nbd[j], dim=1)): + added_ind.append(j) + if len(added_ind) > 0: + nbds = torch.cat([nbds, nbd[added_ind]]) + else: + nbds = torch.cat([nbds, nbd]) + #print(x) + #print("nbds", nbds.size()) + return nbds + + +def _cartesian_neighbors(grouped_x, edge_mat_list): + """ + For given vertices, it returns all neighboring vertices on cartesian product of the graphs given by edge_mat_list + :param grouped_x: 1D Tensor + :param edge_mat_list: + :return: 2d tensor in which each row is 1-hamming distance far from x + """ + neighbor_list = [] + #print(len(edge_mat_list)) + for i in range(len(edge_mat_list)): + nbd_i_elm = edge_mat_list[i][int(grouped_x[i])].nonzero().squeeze(1) + nbd_i = grouped_x.repeat((nbd_i_elm.numel(), 1)) + nbd_i[:, i] = nbd_i_elm + neighbor_list.append(nbd_i) + + return torch.cat(neighbor_list, dim=0) diff --git a/acquisition/acquisition_optimizers/greedy_ascent.py b/acquisition/acquisition_optimizers/greedy_ascent.py new file mode 100755 index 0000000..2737802 --- /dev/null +++ b/acquisition/acquisition_optimizers/greedy_ascent.py @@ -0,0 +1,43 @@ +import torch + +from acquisition.acquisition_functions import expected_improvement +from acquisition.acquisition_marginalization import acquisition_expectation +from acquisition.acquisition_optimizers.graph_utils import neighbors + + +def greedy_ascent(x_init, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func=expected_improvement, max_n_ascent=float('inf'), reference=None): + """ + In order to find local maximum of an acquisition function, at each vertex, + it follows the most increasing edge starting from an initial point + if MAX_ASCENT is infinity, this method tries to find local maximum, otherwise, + it may stop at a noncritical vertex (this option is for a computational reason) + :param x_init: 1d tensor + :param inference_samples: + :param edge_mat_samples: + :param n_vertices: 1D np.array + :param acquisition_func: + :param max_n_ascent: + :param reference: + :return: 1D Tensor, numeric(float) + """ + n_ascent = 0 + x = x_init + #print("x_init shape", x_init.size()) + max_acquisition = acquisition_expectation(x, inference_samples, partition_samples, n_vertices, + acquisition_func, reference) + while n_ascent < max_n_ascent: + x_nbds = neighbors(x[:inference_samples[0].model.kernel.num_discrete], partition_samples, edge_mat_samples, n_vertices, uniquely=True) + x_nbds = torch.cat((x_nbds, x[inference_samples[0].model.kernel.num_discrete:].unsqueeze(0).repeat(x_nbds.size(0), 1)), dim=1) + nbds_acquisition = acquisition_expectation(x_nbds, inference_samples, partition_samples, n_vertices, + acquisition_func, reference) + max_nbd_acquisition, max_nbd_ind = torch.max(nbds_acquisition, 0) + if max_nbd_acquisition > max_acquisition: + max_acquisition = max_nbd_acquisition + x = x_nbds[max_nbd_ind.item()] + n_ascent += 1 + else: + break + #print("x", x) + #print("max", max_acquisition.item()) + return x, max_acquisition.item() diff --git a/acquisition/acquisition_optimizers/simulated_annealing.py b/acquisition/acquisition_optimizers/simulated_annealing.py new file mode 100755 index 0000000..91cde01 --- /dev/null +++ b/acquisition/acquisition_optimizers/simulated_annealing.py @@ -0,0 +1,79 @@ +import math +import random +import time + +import numpy as np + +from simanneal import Annealer +from simanneal.anneal import round_figures + +from acquisition.acquisition_functions import expected_improvement +from acquisition.acquisition_marginalization import acquisition_expectation +from acquisition.acquisition_optimizers.graph_utils import neighbors + + +TMP_FILE_NAME = '' + + +class GraphSimulatedAnnealing(Annealer): + + def __init__(self, initial_state, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func=expected_improvement, reference=None): + """ + + :param initial_state: 1D Tensor + :param inference_samples: + :param partition_samples: + :param edge_mat_samples: + :param n_vertices: + :param acquisition_func: + :param reference: + """ + super(GraphSimulatedAnnealing, self).__init__(initial_state) + self.inference_samples = inference_samples + self.partition_samples = partition_samples + self.edge_mat_samples = edge_mat_samples + self.n_vertices = n_vertices + self.acquisition_func = acquisition_func + self.reference = reference + self.state_history = [] + self.eval_history = [] + + def move(self): + nbds = neighbors(self.state, self.partition_samples, self.edge_mat_samples, self.n_vertices, uniquely=False) + self.state = nbds[np.random.randint(0, nbds.size(0))] + + def energy(self): + # anneal() minimize + evaluation = -acquisition_expectation(self.state, self.inference_samples, self.partition_samples, + self.n_vertices, self.acquisition_func, self.reference).item() + self.state_history.append(self.state.clone()) + self.eval_history.append(evaluation) + return evaluation + + # To overwrite unnecessary printing + def update(self, *args, **kwargs): + pass + + +def simulated_annealing(x_init, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func, reference=None): + """ + Note that Annealer.anneal() MINIMIZES an objective. + :param x_init: + :param inference_samples: + :param partition_samples: + :param edge_mat_samples: + :param n_vertices: + :param acquisition_func: + :param reference: + :return: 1D Tensor, numeric(float) + """ + sa_runner = GraphSimulatedAnnealing(x_init, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func, reference) + steps = 500 + sa_runner.set_schedule({'tmax': 1.0, 'tmin': 0.8 ** steps, 'steps': steps, 'updates': sa_runner.updates}) + opt_state, opt_eval = sa_runner.anneal() + + # Annealer.anneal() MINinimzes an objective but acqusition functions should be MAXimized. + return opt_state, -opt_eval diff --git a/acquisition/acquisition_optimizers/starting_points.py b/acquisition/acquisition_optimizers/starting_points.py new file mode 100755 index 0000000..145211e --- /dev/null +++ b/acquisition/acquisition_optimizers/starting_points.py @@ -0,0 +1,62 @@ +import numpy as np + +import torch + +from acquisition.acquisition_optimizers.graph_utils import neighbors +from acquisition.acquisition_marginalization import acquisition_expectation +from acquisition.acquisition_functions import expected_improvement + + +N_RANDOM_VERTICES = 200#20000 +N_GREEDY_ASCENT_INIT = 20 +N_SPRAY = 10 + + +def optim_inits(objective, x_opt, inference_samples, partition_samples, edge_mat_samples, n_vertices, + acquisition_func=expected_improvement, reference=None): + """ + :param x_opt: 1D Tensor + :param inference_samples: + :param partition_samples: + :param edge_mat_samples: + :param n_vertices: + :param acquisition_func: + :param reference: + :return: + """ + # for x, y in zip(objective.problem.lower_bounds, objective.problem.upper_bounds): + # print(x, y) + # print("n_vertices", n_vertices) + # print(partition_samples) + #print(edge_mat_samples) + #print(len(edge_mat_samples)) + #print(edge_mat_samples[0]) + #for i in range(len(edge_mat_samples)): + # print(len(edge_mat_samples[i])) + + #rnd_nbd = torch.cat(tuple([torch.randint(low=0, high=int(n_v), size=(N_RANDOM_VERTICES, 1)) for n_v in n_vertices]), dim=1).long() + rnd_nbd = objective.generate_random_points(N_RANDOM_VERTICES) + min_nbd = neighbors(x_opt[:objective.num_discrete], partition_samples, edge_mat_samples, n_vertices, uniquely=False) + # print(min_nbd.size(0)) + # print(min_nbd) + # print(x_opt[objective.num_discrete:].unsqueeze(0).repeat(min_nbd.size(0), 1)[:10]) + min_nbd = torch.cat((min_nbd, x_opt[objective.num_discrete:].unsqueeze(0).repeat(min_nbd.size(0), 1)), dim=1) + # print(min_nbd[:6]) + shuffled_ind = list(range(min_nbd.size(0))) + np.random.shuffle(shuffled_ind) + x_init_candidates = torch.cat(tuple([min_nbd[shuffled_ind[:N_SPRAY]], rnd_nbd]), dim=0) + acquisition_values = acquisition_expectation(x_init_candidates, inference_samples, partition_samples, n_vertices, + acquisition_func, reference) + #print("acquisition_values") + #print(acquisition_values[:30]) + + nonnan_ind = ~torch.isnan(acquisition_values).squeeze(1) + x_init_candidates = x_init_candidates[nonnan_ind] + acquisition_values = acquisition_values[nonnan_ind] + + acquisition_sorted, acquisition_sort_ind = torch.sort(acquisition_values.squeeze(1), descending=True) + x_init_candidates = x_init_candidates[acquisition_sort_ind] + + return x_init_candidates[:N_GREEDY_ASCENT_INIT], acquisition_sorted[:N_GREEDY_ASCENT_INIT] + + diff --git a/config.py b/config.py new file mode 100755 index 0000000..03ded59 --- /dev/null +++ b/config.py @@ -0,0 +1,8 @@ +PROGRESS_BAR_LEN = 50 + +def experiment_directory(): + return '../EXPERIMENTS' + + +def data_directory(): + return '../EXPERIMENTS' diff --git a/experiments/__init__.py b/experiments/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/experiments/exp_utils.py b/experiments/exp_utils.py new file mode 100755 index 0000000..cd46ad5 --- /dev/null +++ b/experiments/exp_utils.py @@ -0,0 +1,25 @@ +import torch + + +def sample_init_points(n_vertices, n_points, random_seed=None): + """ + + :param n_vertices: 1D np.array + :param n_points: + :param random_seed: + :return: + """ + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = torch.empty(0).long() + for _ in range(n_points): + init_points = torch.cat([init_points, torch.cat([torch.randint(0, int(elm), (1, 1)) for elm in n_vertices], dim=1)], dim=0) + if random_seed is not None: + torch.set_rng_state(rng_state) + return init_points + + +if __name__ == '__main__': + print(sample_init_points([2] * 10, 5, 3)) + diff --git a/experiments/random_seed_config.py b/experiments/random_seed_config.py new file mode 100755 index 0000000..bde2d25 --- /dev/null +++ b/experiments/random_seed_config.py @@ -0,0 +1,22 @@ +import numpy as np + + +SEED_STR_LIST = ['2021AAAI_COCO'] + + +def generate_random_seed_coco(): + return _generate_random_seed('2021AAAI_COCO', n_init_point_seed=25) + +def _generate_random_seed(seed_str, n_init_point_seed=10): + assert seed_str in SEED_STR_LIST + rng_state = np.random.RandomState(seed=sum([ord(ch) for ch in seed_str])) + return rng_state.randint(0, 10000, (n_init_point_seed, )) + + +def _generate_random_seed_pair(seed_str, n_test_case_seed=5, n_init_point_seed=5): + assert seed_str in SEED_STR_LIST + rng_state = np.random.RandomState(seed=sum([ord(ch) for ch in seed_str])) + result = {} + for _ in range(n_test_case_seed): + result[rng_state.randint(0, 10000)] = list(rng_state.randint(0, 10000, (n_init_point_seed, ))) + return result diff --git a/experiments/test_functions/__init__.py b/experiments/test_functions/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/experiments/test_functions/em_func.py b/experiments/test_functions/em_func.py new file mode 100644 index 0000000..26baffe --- /dev/null +++ b/experiments/test_functions/em_func.py @@ -0,0 +1,100 @@ +from __future__ import division, print_function +import numpy as np +import torch +class Problem(object): + def __init__(self, dimension, lower_bounds, upper_bounds): + self.dimension = dimension + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + pass + def c(self, s, t, M, D, L, tau): + print((M*np.exp(-(s**2)/4*D*t))/np.sqrt(4*np.pi*D*t) + ((t > tau) * M * np.exp(-((s-L)**2)/(4*D*(t-tau))))/np.sqrt(4*np.pi*D*(t-tau))) + val = (M*np.exp(-(s**2)/4*D*t))/np.sqrt(4*np.pi*D*t) + if (t > tau): + val += ((t > tau) * M * np.exp(-((s-L)**2)/(4*D*(t-tau))))/np.sqrt(4*np.pi*D*(t-tau)) + return val + def __call__(self, x): + tau = 30.01 + x[0]/1000 + M = 7 + ((x[1] + 1)*6)/2 + D = 0.02 + ((x[2] + 1)*0.10)/2 + L = 0.01 + ((x[3] + 1)*2.99)/2 + print(f"M:{M}, D:{D}, L:{L}, tau:{tau}") + val = 0.0 + for s in [0, 1, 2.5]: + for t in [15, 30, 45, 60]: + print(f"s:{s}, t:{t}") + val += (self.c(s, t, 10, 0.07, 1.505, 30.1525) - self.c(s, t, M, D, L, tau))**2 + print(f"val:{val}") + return val +class EM_func(object): + """ + Environmental model function + """ + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + self.lower_bounds = [] + self.upper_bounds = [] + # discrete variables + self.lower_bounds.append(0) + self.upper_bounds.append(284) + for i in range(3): + self.lower_bounds.append(-1) + self.upper_bounds.append(1) + assert len(self.lower_bounds) == 4 + assert len(self.upper_bounds) == 4 + #print(lower_bounds) + #print(upper_bounds) + self.num_discrete = 1 + self.num_continuous = 3 + + self.problem = Problem(dimension=4, lower_bounds=self.lower_bounds, upper_bounds=self.upper_bounds) + self.problem.dimension = 4 + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [] + for i in range(self.num_discrete): + self.n_vertices.append(self.upper_bounds[i] - self.lower_bounds[i] + 1) + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.generate_random_points(n_points=10, random_seed=random_seed).float() + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x= x_unorder.numpy().copy() + print(f"evaluating {x}....") + evaluation = self.problem(x) + print(evaluation) + return torch.tensor(evaluation).float() + + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + for i in range(self.num_discrete): + random_point.append(torch.randint(self.lower_bounds[i], self.upper_bounds[i]+1, (1,))) + for i in range(self.num_discrete, self.num_discrete + self.num_continuous): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + diff --git a/experiments/test_functions/experiment_configuration.py b/experiments/test_functions/experiment_configuration.py new file mode 100755 index 0000000..aac909a --- /dev/null +++ b/experiments/test_functions/experiment_configuration.py @@ -0,0 +1,155 @@ +import os +import numpy as np +import scipy.io as sio + +import torch + + +ISING_GRID_H = 4 +ISING_GRID_W = 4 +ISING_N_EDGES = 24 + +CONTAMINATION_N_STAGES = 25 + +AEROSTRUCTURAL_N_COUPLINGS = 21 + + +def sample_speed_reducer_points(n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + random_point.append(torch.randint(0, 12, (1,))) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + random_point.append(torch.FloatTensor(1).uniform_(0, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + +def sample_weld_points(n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + random_point.append(torch.randint(0, 2, (1,))) + random_point.append(torch.randint(0, 4, (1,))) + random_point.append(torch.FloatTensor(1).uniform_(0.0625, 2)) + random_point.append(torch.FloatTensor(1).uniform_(0, 10)) + random_point.append(torch.FloatTensor(1).uniform_(2, 10)) + random_point.append(torch.FloatTensor(1).uniform_(0.0625, 2)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + +def sample_mixed_init_points(lbs, ubs, num_discrete, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + for i in range(num_discrete): + random_point.append(torch.randint(int(lbs[i]), int(ubs[i]+1), (1,))) + for i in range(num_discrete, len(lbs)): + random_point.append(torch.FloatTensor(1).uniform_(lbs[i], ubs[i])) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + + +def sample_init_points(n_vertices, n_points, random_seed=None): + """ + + :param n_vertices: 1D array + :param n_points: + :param random_seed: + :return: + """ + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = torch.empty(0).long() + for _ in range(n_points): + init_points = torch.cat([init_points, torch.cat([torch.randint(0, int(elm), (1, 1)) for elm in n_vertices], dim=1)], dim=0) + if random_seed is not None: + torch.set_rng_state(rng_state) + return init_points + + +def generate_ising_interaction(grid_h, grid_w, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + horizontal_interaction = ((torch.randint(0, 2, (grid_h * (grid_w - 1), )) * 2 - 1).float() * (torch.rand(grid_h * (grid_w - 1)) * (5 - 0.05) + 0.05)).view(grid_h, grid_w-1) + vertical_interaction = ((torch.randint(0, 2, ((grid_h - 1) * grid_w, )) * 2 - 1).float() * (torch.rand((grid_h - 1) * grid_w) * (5 - 0.05) + 0.05)).view(grid_h-1, grid_w) + if random_seed is not None: + torch.set_rng_state(rng_state) + return horizontal_interaction, vertical_interaction + + +def generate_contamination_dynamics(random_seed=None): + n_stages = CONTAMINATION_N_STAGES + n_simulations = 100 + + init_alpha = 1.0 + init_beta = 30.0 + contam_alpha = 1.0 + contam_beta = 17.0 / 3.0 + restore_alpha = 1.0 + restore_beta = 3.0 / 7.0 + init_Z = np.random.RandomState(random_seed).beta(init_alpha, init_beta, size=(n_simulations,)) + lambdas = np.random.RandomState(random_seed).beta(contam_alpha, contam_beta, size=(n_stages, n_simulations)) + gammas = np.random.RandomState(random_seed).beta(restore_alpha, restore_beta, size=(n_stages, n_simulations)) + + return init_Z, lambdas, gammas + + +def interaction_sparse2dense(bocs_representation): + assert bocs_representation.size(0) == bocs_representation.size(1) + grid_size = int(bocs_representation.size(0) ** 0.5) + horizontal_interaction = torch.zeros(grid_size, grid_size-1) + vertical_interaction = torch.zeros(grid_size-1, grid_size) + for i in range(bocs_representation.size(0)): + r_i = i // grid_size + c_i = i % grid_size + for j in range(i + 1, bocs_representation.size(1)): + r_j = j // grid_size + c_j = j % grid_size + if abs(r_i - r_j) + abs(c_i - c_j) > 1: + assert bocs_representation[i, j] == 0 + elif abs(r_i - r_j) == 1: + vertical_interaction[min(r_i, r_j), c_i] = bocs_representation[i, j] + else: + horizontal_interaction[r_i, min(c_i, c_j)] = bocs_representation[i, j] + return horizontal_interaction, vertical_interaction + + +def interaction_dense2sparse(horizontal_interaction, vertical_interaction): + grid_size = horizontal_interaction.size(0) + bocs_representation = torch.zeros(grid_size ** 2, grid_size ** 2) + for i in range(bocs_representation.size(0)): + r_i = i // grid_size + c_i = i % grid_size + for j in range(i + 1, bocs_representation.size(1)): + r_j = j // grid_size + c_j = j % grid_size + if abs(r_i - r_j) + abs(c_i - c_j) > 1: + assert bocs_representation[i, j] == 0 + elif abs(r_i - r_j) == 1: + bocs_representation[i, j] = vertical_interaction[min(r_i, r_j), c_i] + else: + bocs_representation[i, j] = horizontal_interaction[r_i, min(c_i, c_j)] + return bocs_representation + bocs_representation.t() + + +if __name__ == '__main__': + _convert_random_data_to_matfile() diff --git a/experiments/test_functions/mixed_integer.py b/experiments/test_functions/mixed_integer.py new file mode 100644 index 0000000..38cd346 --- /dev/null +++ b/experiments/test_functions/mixed_integer.py @@ -0,0 +1,62 @@ +from __future__ import division, print_function +from experiments.test_functions.experiment_configuration import sample_mixed_init_points +import cocoex, cocopp # experimentation and post-processing modules + +# prepare mixed integer suite +suite_name = "bbob-mixint" +output_folder = "cocex-optimize-fmin" +suite = cocoex.Suite(suite_name, "", "") +#observer = cocoex.Observer(suite_name, "result_folder: " + output_folder) +#minimal_print = cocoex.utilities.MiniPrint() +import numpy as np +import torch + +class MixedIntegerCOCO(object): + """ + Mixed Integer Black box optimization using cocoex library + """ + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + self.problem = suite.get_problem(self.problem_id) + self.num_discrete = self.problem.number_of_integer_variables + self.num_continuous = self.problem.dimension - self.num_discrete + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [] + for i in range(self.num_discrete): + self.n_vertices.append(int(self.problem.upper_bounds[i] - self.problem.lower_bounds[i] + 1)) + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = sample_mixed_init_points(self.problem.lower_bounds, self.problem.upper_bounds, + self.num_discrete, n_points=20, random_seed=random_seed).float() + self.suggested_init[-1] = torch.tensor(self.problem.initial_solution).float() + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x): + if x.dim() == 2: + x = x.squeeze(0) + evaluation = self.problem(x) + return torch.tensor(evaluation).float() + + def generate_random_points(self, n_points, random_seed=None): + return sample_mixed_init_points(self.problem.lower_bounds, self.problem.upper_bounds, + self.num_discrete, n_points=n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + +if __name__ == '__main__': + mixobj = MixedIntegerCOCO(44,'bbob-mixint_f001_i01_d20') + print(mixobj.evaluate(mixobj.suggested_init[0])) + print(mixobj.problem.evaluations) + print(mixobj.evaluate(mixobj.suggested_init[5])) + print(mixobj.problem.evaluations) diff --git a/experiments/test_functions/ml_datasets.py b/experiments/test_functions/ml_datasets.py new file mode 100644 index 0000000..9d41687 --- /dev/null +++ b/experiments/test_functions/ml_datasets.py @@ -0,0 +1,87 @@ +import numpy as np +import timeit +from sklearn import datasets +from sklearn.datasets import fetch_olivetti_faces +import openml +from sklearn.model_selection import train_test_split, cross_val_score +from sklearn.preprocessing import MinMaxScaler +from sklearn.naive_bayes import BernoulliNB, MultinomialNB +from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis +from sklearn.svm import LinearSVC, SVC +from sklearn.linear_model import SGDClassifier, LogisticRegression +from sklearn.linear_model.passive_aggressive import PassiveAggressiveClassifier +from sklearn.neural_network import MLPClassifier +from sklearn.tree import DecisionTreeClassifier +from sklearn.ensemble import AdaBoostClassifier, ExtraTreesClassifier, GradientBoostingClassifier, RandomForestClassifier +from sklearn.metrics import accuracy_score + +datasets_all = ["iris", "digits", "wine", "breast_cancer", "analcatdata_authorship", + "blood_transfusion", "monks1", "monks2", "steel_plates_fault", "qsar_biodeg", + "phoneme", "diabetes", "hill_valley", "eeg_eye_state", "waveform", + "spambase", "australian", "churn", "vehicle", "balance_scale", + "kc1", "kc2", "cardiotocography", "wall_robot_navigation", "segment", + "artificial_characters", "electricity", "gas_drift", "olivetti", "letter"] + +# cross-validation on training data +cv_train = 5 + +def gen_train_test_data(dataset="", seed=42): + print("dataset: {}, seed: {}".format(dataset, seed)) + if dataset == "blood_transfusion": # 748 + dataset_id = 1464 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + if dataset == "phoneme": # 5404 + dataset_id = 1489 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + if dataset == "kc1": # 2109 + dataset_id = 1067 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + y_new = np.array([str(val) for val in y.values]) + y_new[y_new == "False"] = "0" + y_new[y_new == "True"] = "1" + y = np.array([int(val) for val in y_new]) + if dataset == "australian": # 4202 + dataset_id = 40981 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + if dataset == "vehicle": # 846 + dataset_id = 54 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + y_new = np.array([val for val in y.values]) + y_new[y_new == "opel"] = "0" + y_new[y_new == "saab"] = "1" + y_new[y_new == "bus"] = "2" + y_new[y_new == "van"] = "3" + y = np.array([int(val) for val in y_new]) + if dataset == "connect4": + dataset_id = 40668 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + if dataset == "segment": + dataset_id = 40984 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + y_new = np.array([val for val in y.values]) + y_new[y_new == "brickface"] = "0" + y_new[y_new == "sky"] = "1" + y_new[y_new == "foliage"] = "2" + y_new[y_new == "cement"] = "3" + y_new[y_new == "window"] = "4" + y_new[y_new == "path"] = "5" + y_new[y_new == "grass"] = "6" + y = np.array([int(val) for val in y_new]) + if dataset == "cnae": + dataset_id = 1468 + ds = openml.datasets.get_dataset(dataset_id) + (X, y, categorical, names) = ds.get_data(target=ds.default_target_attribute) + # normalize X + X = MinMaxScaler().fit_transform(X) + # convert y to integer array + y = np.array([int(val) for val in y]) + X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=seed) + print(f"X_train shape: {X_train.shape}, X_test shape: {X_test.shape}") + return X_train, y_train, X_test, y_test diff --git a/experiments/test_functions/nn_ml_datasets.py b/experiments/test_functions/nn_ml_datasets.py new file mode 100644 index 0000000..abf9294 --- /dev/null +++ b/experiments/test_functions/nn_ml_datasets.py @@ -0,0 +1,126 @@ +from __future__ import division, print_function +import numpy as np +import torch +from sklearn.neural_network import MLPClassifier +from sklearn.metrics import mean_squared_error, accuracy_score +import experiments.test_functions.ml_datasets as mld + +class Problem(object): + def __init__(self, dataset, dimension, num_discrete, num_continuous, lower_bounds, upper_bounds): + self.dataset = dataset + self.dimension = dimension + self.num_discrete = num_discrete + self.num_continuous = num_continuous + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + self.train_X, self.train_y, self.test_X, self.test_y = mld.gen_train_test_data(self.dataset) + + def __call__(self, x): + hidden_layer_sizes = list(range(40, 310, 20)) + activation = ['identity', 'logistic', 'tanh', 'relu'] + batch_size = list(range(40, 210, 20)) + learning_rate = ['constant', 'invscaling', 'adaptive'] + early_stopping = [False, True] + learning_rate_init = [0.001, 1] + momentum = [0.5, 1] + alpha = [0.0001, 1] + x[5] = 0.001 + ((x[5] + 1)*(1-0.001))/2 + x[6] = 0.5 + ((x[6] + 1)*(1-0.5))/2 + x[7] = 0.0001 + ((x[7] + 1)*(1-0.0001))/2 + print(f"x after conversion {x}") + kwag_dict = {'hidden_layer_sizes': hidden_layer_sizes[int(x[0])], 'activation': activation[int(x[1])], 'batch_size': batch_size[int(x[2])], 'learning_rate': learning_rate[int(x[3])], + 'early_stopping': early_stopping[int(x[4])], 'learning_rate_init': x[5], 'momentum':x[6], 'alpha':x[7], 'solver':'sgd'} + print(kwag_dict) + model = MLPClassifier(**kwag_dict, random_state=2) + model.fit(self.train_X, self.train_y) + pred_y = model.predict(self.test_X) + # maximize accuracy + auc = accuracy_score(self.test_y, pred_y) + print(f"auc: {auc}") + return -auc + +class NN_ML_Datasets(object): + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + self.lower_bounds = [] + self.upper_bounds = [] + # Discrete bounds + self.lower_bounds.append(0) + self.upper_bounds.append(13) + self.lower_bounds.append(0) + self.upper_bounds.append(3) + + self.lower_bounds.append(0) + self.upper_bounds.append(8) + + self.lower_bounds.append(0) + self.upper_bounds.append(2) + self.lower_bounds.append(0) + self.upper_bounds.append(1) + + + # Continuous bounds + for i in range(3): + self.lower_bounds.append(-1) + self.upper_bounds.append(1) + assert len(self.lower_bounds) == 8 + assert len(self.upper_bounds) == 8 + print(self.lower_bounds) + print(self.upper_bounds) + self.num_discrete = 5 + self.num_continuous = 3 + self.dimension = 8 + assert self.dimension == self.num_discrete + self.num_continuous + self.problem = Problem(dataset=self.problem_id, dimension=self.dimension, num_discrete=self.num_discrete, num_continuous=self.num_continuous, lower_bounds=self.lower_bounds, upper_bounds=self.upper_bounds) + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [] + for i in range(self.num_discrete): + self.n_vertices.append(self.upper_bounds[i] - self.lower_bounds[i] + 1) + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.sample_points(n_points=10, random_seed=random_seed).float() + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + #assert x.numel() == len(self.n_vertices) + #x = x.numpy() + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x= x_unorder.numpy().copy() + print(f"evaluating {x}....") + evaluation = self.problem(x) + #print(evaluation) + #print(type(evaluation)) + return torch.tensor(evaluation).float() + + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + for i in range(self.num_discrete): + random_point.append(torch.randint(self.lower_bounds[i], self.upper_bounds[i] + 1, (1,))) + for i in range(self.num_discrete, self.num_discrete+self.num_continuous): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + diff --git a/experiments/test_functions/pressure_vessel_design.py b/experiments/test_functions/pressure_vessel_design.py new file mode 100644 index 0000000..b98e7ec --- /dev/null +++ b/experiments/test_functions/pressure_vessel_design.py @@ -0,0 +1,93 @@ +from __future__ import division, print_function +import numpy as np +import torch +class Problem(object): + def __init__(self, dimension, lower_bounds, upper_bounds): + self.dimension = dimension + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + pass + def __call__(self, x): + x[0] += 1 + x[1] += 1 + x[2] = 10 + ((x[2] + 1) * (200 - 10))/2 + x[3] = 10 + ((x[3] + 1) * (240 - 10))/2 + func_value = 0.6224 * x[0] * x[2] * x[3] + 1.7781 * x[1] * (x[2] ** 2) + 3.1661 * (x[0] ** 2) * x[3] + 19.84 * (x[0] ** 2) * x[2] + constraint_value = 0 + constraint_1 = x[0] - 0.0193*x[2] + constraint_2 = x[1] - 0.00954*x[2] + constraint_3 = np.pi*(x[2]**2)*x[3] + (4/3)*np.pi*(x[2] ** 3) - 1296000 + constraint_cost = (constraint_1 < 0) * 100 + (constraint_2 < 0) * 100 + (constraint_3 < 0) * 100 + print(f"func_value : {func_value/1e6}") + return func_value/1e6 #+ constraint_cost + +class Pressure_Vessel_Design(object): + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + lower_bounds = [] + upper_bounds = [] + # discrete variables + lower_bounds.append(0) # + upper_bounds.append(99) + lower_bounds.append(0) # + upper_bounds.append(99) + + for i in range(2): + lower_bounds.append(-1) # weld_thickness (h) + upper_bounds.append(1) + + assert len(lower_bounds) == 4 + assert len(upper_bounds) == 4 + print(lower_bounds) + print(upper_bounds) + self.problem = Problem(dimension=4, lower_bounds=lower_bounds, upper_bounds=upper_bounds) + self.problem.dimension = 4 + self.num_discrete = 2 + self.num_continuous = 2 + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [100, 100] + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.generate_random_points(n_points=10, random_seed=random_seed) + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x = x_unorder.numpy().copy() + print(f"evaluating {x}....") + evaluation = self.problem(x) + return torch.tensor(evaluation).float() + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + random_point.append(torch.randint(0, 100, (1,))) + random_point.append(torch.randint(0, 100, (1,))) + + for i in range(2): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + + diff --git a/experiments/test_functions/push_robot_14d.py b/experiments/test_functions/push_robot_14d.py new file mode 100644 index 0000000..7e13d4b --- /dev/null +++ b/experiments/test_functions/push_robot_14d.py @@ -0,0 +1,108 @@ +from __future__ import division, print_function +import numpy as np +import torch +from experiments.test_functions.robot_push_14d import push_function +from hyperopt import hp +class Problem(object): + def __init__(self, dimension, lower_bounds, upper_bounds): + self.dimension = dimension + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + pass + def __call__(self, x): + argv = [] + argv.append(x[0]-5) + argv.append(x[1]-5) + argv.append(x[4]-10) + argv.append(x[5]-10) + argv.append(x[8]+2) + argv.append(((x[10]+1)*2*np.pi)/2) + argv.append(x[2]-5) + argv.append(x[3]-5) + argv.append(x[6]-10) + argv.append(x[7]-10) + argv.append(x[9]+2) + argv.append(((x[11]+1)*2*np.pi)/2) + argv.append(((x[12]+1)*10)/2) + argv.append(((x[13]+1)*10)/2) + f = push_function.PushReward() + reward = f(argv) + print(f"reward: {reward}") + return -1*reward + +class Push_robot_14d(object): + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + self.lower_bounds = [] + self.upper_bounds = [] + # discrete variables + for i in range(4): + self.lower_bounds.append(0) + self.upper_bounds.append(10) + for i in range(4): + self.lower_bounds.append(0) + self.upper_bounds.append(20) + for i in range(2): + self.lower_bounds.append(0) + self.upper_bounds.append(28) + for i in range(4): + self.lower_bounds.append(-1) + self.upper_bounds.append(1) + assert len(self.lower_bounds) == 14 + assert len(self.upper_bounds) == 14 + #print(lower_bounds) + #print(upper_bounds) + self.num_discrete = 10 + self.num_continuous = 4 + + self.problem = Problem(dimension=14, lower_bounds=self.lower_bounds, upper_bounds=self.upper_bounds) + self.problem.dimension = 14 + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [] + for i in range(self.num_discrete): + self.n_vertices.append(self.upper_bounds[i] - self.lower_bounds[i] + 1) + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.generate_random_points(n_points=10, random_seed=random_seed).float() + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x= x_unorder.numpy().copy() + print(f"evaluating {x}....") + evaluation = self.problem(x) + print(evaluation) + return torch.tensor(evaluation).float() + + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + for i in range(self.num_discrete): + random_point.append(torch.randint(self.lower_bounds[i], self.upper_bounds[i], (1,))) + for i in range(self.num_discrete, self.num_discrete + self.num_continuous): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + diff --git a/experiments/test_functions/push_world.py b/experiments/test_functions/push_world.py new file mode 100644 index 0000000..8e137e9 --- /dev/null +++ b/experiments/test_functions/push_world.py @@ -0,0 +1,275 @@ +#!/usr/bin/env python +from Box2D import * +from Box2D.b2 import * +import numpy as np +import pygame +import scipy.io +from numpy import linalg as LA + +# this just makes pygame show what's going on +class guiWorld: + def __init__(self, fps): + self.SCREEN_WIDTH, self.SCREEN_HEIGHT = 1000, 1000 + self.TARGET_FPS = fps + self.PPM = 10.0 # pixels per meter + self.screen = pygame.display.set_mode((self.SCREEN_WIDTH, self.SCREEN_HEIGHT), 0, 32) + pygame.display.set_caption('push simulator') + self.clock = pygame.time.Clock() + self.screen_origin = b2Vec2(self.SCREEN_WIDTH/(2*self.PPM), self.SCREEN_HEIGHT/(self.PPM*2)) + self.colors = { + b2_staticBody : (255,255,255,255), + b2_dynamicBody : (163,209,224,255) + } + + def draw(self, bodies, bg_color=(64,64,64,0)): + #def draw(self, bodies, bg_color=(0,0,0,0)): + def my_draw_polygon(polygon, body, fixture): + vertices=[(self.screen_origin + body.transform*v)*self.PPM for v in polygon.vertices] + vertices=[(v[0], self.SCREEN_HEIGHT-v[1]) for v in vertices] + color = self.colors[body.type] + if body.userData == "obs": + color = (123,128,120,0) + if body.userData == "hand": + color = (174,136,218,0) + + pygame.draw.polygon(self.screen, color, vertices) + + def my_draw_circle(circle, body, fixture): + position=(self.screen_origin + body.transform*circle.pos)*self.PPM + position=(position[0], self.SCREEN_HEIGHT-position[1]) + color = self.colors[body.type] + if body.userData == "hand": + color = (174,136,218,0) + pygame.draw.circle(self.screen, color, [int(x) for x in + position], int(circle.radius*self.PPM)) + + b2PolygonShape.draw=my_draw_polygon + b2CircleShape.draw=my_draw_circle + # draw the world + self.screen.fill(bg_color) + self.clock.tick(self.TARGET_FPS) + for body in bodies: + for fixture in body.fixtures: + fixture.shape.draw(body,fixture) + pygame.display.flip() + +# this is the interface to pybox2d +class b2WorldInterface: + def __init__(self, do_gui=True): + self.world = b2World(gravity=(0.0,0.0), doSleep=True) + self.do_gui = do_gui + self.TARGET_FPS = 100 + self.TIME_STEP = 1.0/self.TARGET_FPS + self.VEL_ITERS, self.POS_ITERS =10,10 + self.bodies = [] + + if do_gui: + self.gui_world = guiWorld(self.TARGET_FPS) + #raw_input() + else: + self.gui_world = None + + def initialize_gui(self): + if self.gui_world == None: + self.gui_world = guiWorld(self.TARGET_FPS) + self.do_gui = True + def stop_gui(self): + self.do_gui = False + + def add_bodies(self, new_bodies): + """ add a single b2Body or list of b2Bodies to the world""" + if type(new_bodies) == list: + self.bodies += new_bodies + else: + self.bodies.append(new_bodies) + def step(self, show_display=True, idx=0): + self.world.Step(self.TIME_STEP, self.VEL_ITERS, self.POS_ITERS) + if show_display and self.do_gui: + self.gui_world.draw(self.bodies) + #if idx % 10 == 0: + # pygame.image.save(self.gui_world.screen,'tmp_images/'+str(int(sm.ttt*100)+idx)+'.bmp') + +class end_effector: + def __init__(self, b2world_interface, init_pos, base, init_angle, hand_shape='rectangle', hand_size=(0.3,1)): + world= b2world_interface.world + self.hand = world.CreateDynamicBody(position=init_pos,angle=init_angle) + self.hand_shape = hand_shape + self.hand_size = hand_size + # forceunit for circle and rect + if hand_shape == 'rectangle': + rshape = b2PolygonShape(box=hand_size) + self.forceunit = 30.0 + elif hand_shape == 'circle': + rshape = b2CircleShape(radius=hand_size) + self.forceunit = 100.0 + elif hand_shape == 'polygon': + rshape = b2PolygonShape(vertices=hand_size) + else: + raise Exception("%s is not a correct shape" % hand_shape) + + self.hand.CreateFixture( + shape = rshape, + density = .1, + friction = .1 + ) + self.hand.userData = "hand" + + friction_joint = world.CreateFrictionJoint( + bodyA = base, + bodyB = self.hand, + maxForce = 2, + maxTorque = 2, + ) + b2world_interface.add_bodies(self.hand) + + + + def set_pos(self, pos, angle): + self.hand.position = pos + self.hand.angle = angle + def apply_wrench(self, rlvel=(0,0), ravel=0): + #self.hand.ApplyForce(force, self.hand.position,wake=True) + #if avel != 0: + + avel = self.hand.angularVelocity + delta_avel = ravel - avel + torque = self.hand.mass*delta_avel*30.0 + self.hand.ApplyTorque(torque, wake=True) + + #else: + lvel = self.hand.linearVelocity + delta_lvel = b2Vec2(rlvel) - b2Vec2(lvel) + force = self.hand.mass*delta_lvel*self.forceunit + self.hand.ApplyForce(force, self.hand.position,wake=True) + + + def get_state(self, verbose=False): + state = list(self.hand.position) + [ self.hand.angle] + \ + list(self.hand.linearVelocity) + [self.hand.angularVelocity] + if verbose: + print_state = ["%.3f" % x for x in state] + #print "position, velocity: (%s), (%s) " % \ + # ((", ").join(print_state[:3]), (", ").join(print_state[3:]) ) + + return state + +def make_thing(table_width, table_length, b2world_interface, thing_shape, thing_size, thing_friction, thing_density, base_friction, obj_loc): + world = b2world_interface.world + base = world.CreateStaticBody( + position = (0,0), + #friction = base_friction, + shapes = b2PolygonShape(box=(table_length,table_width)), + ) + + link = world.CreateDynamicBody(position=obj_loc) + if thing_shape == 'rectangle': + linkshape = b2PolygonShape(box=thing_size) + elif thing_shape == 'circle': + linkshape = b2CircleShape(radius=thing_size) + elif thing_shape == 'polygon': + linkshape = b2PolygonShape(vertices=thing_size) + else: + raise Exception("%s is not a correct shape" % thing_shape) + + link.CreateFixture( + shape = linkshape, + density = thing_density, + friction = thing_friction, + ) + friction_joint = world.CreateFrictionJoint( + bodyA = base, + bodyB = link, + maxForce = 5, + maxTorque = 2, + ) + + b2world_interface.add_bodies([base,link]) + return link,base + + +def simu_push(world, thing, robot, base, simulation_steps): + # simulating push with fixed direction pointing from robot location to thing location + desired_vel = thing.position - robot.hand.position + desired_vel = desired_vel / np.linalg.norm(desired_vel) * 5 + rvel = b2Vec2(desired_vel[0]+np.random.normal(0,0.1),desired_vel[1]+np.random.normal(0,0.1)) + + rstop = False + for t in range(simulation_steps+100): + if not rstop: + robot.apply_wrench(rvel) + world.step() + + ostate = list(thing.position) + [ thing.angle] + \ + list(thing.linearVelocity) + [thing.angularVelocity] + if t == simulation_steps - 1: + rstop = True + return list(thing.position) +def simu_push2(world, thing, robot, base, xvel, yvel, simulation_steps): + desired_vel = np.array([xvel, yvel]) + rvel = b2Vec2(desired_vel[0]+np.random.normal(0,0.1),desired_vel[1]+np.random.normal(0,0.1)) + + rstop = False + for t in range(simulation_steps+100): + if not rstop: + robot.apply_wrench(rvel) + world.step() + + ostate = list(thing.position) + [ thing.angle] + \ + list(thing.linearVelocity) + [thing.angularVelocity] + if t == simulation_steps - 1: + rstop = True + return list(thing.position) + +def make_1thing(base, b2world_interface, thing_shape, thing_size, thing_friction, thing_density, obj_loc): + world = b2world_interface.world + + link = world.CreateDynamicBody(position=obj_loc) + if thing_shape == 'rectangle': + linkshape = b2PolygonShape(box=thing_size) + elif thing_shape == 'circle': + linkshape = b2CircleShape(radius=thing_size) + elif thing_shape == 'polygon': + linkshape = b2PolygonShape(vertices=thing_size) + else: + raise Exception("%s is not a correct shape" % thing_shape) + + link.CreateFixture( + shape = linkshape, + density = thing_density, + friction = thing_friction, + ) + friction_joint = world.CreateFrictionJoint( + bodyA = base, + bodyB = link, + maxForce = 5, + maxTorque = 2, + ) + + b2world_interface.add_bodies([link]) + return link +def simu_push_2robot2thing(world, thing, thing2, robot, robot2, base, xvel, yvel, xvel2, yvel2, rtor, rtor2, simulation_steps, simulation_steps2): + desired_vel = np.array([xvel, yvel]) + rvel = b2Vec2(desired_vel[0]+np.random.normal(0,1e-6),desired_vel[1]+np.random.normal(0,1e-6)) + + desired_vel2 = np.array([xvel2, yvel2]) + rvel2 = b2Vec2(desired_vel2[0]+np.random.normal(0,1e-6),desired_vel2[1]+np.random.normal(0,1e-6)) + tmax = np.max([simulation_steps,simulation_steps2]) + for t in range(tmax+100): + if t < simulation_steps: + robot.apply_wrench(rvel, rtor) + if t < simulation_steps2: + robot2.apply_wrench(rvel2, rtor2) + world.step() + + return (list(thing.position), list(thing2.position)) +def make_base(table_width, table_length, b2world_interface): + world = b2world_interface.world + base = world.CreateStaticBody( + position = (0,0), + #friction = base_friction, + shapes = b2PolygonShape(box=(table_length,table_width)), + ) + + + b2world_interface.add_bodies([base]) + return base diff --git a/experiments/test_functions/robot_push_14d/__init__.py b/experiments/test_functions/robot_push_14d/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/experiments/test_functions/robot_push_14d/__pycache__/__init__.cpython-36.pyc b/experiments/test_functions/robot_push_14d/__pycache__/__init__.cpython-36.pyc new file mode 100644 index 0000000000000000000000000000000000000000..088f8d77eff85ab50080689cf81aa024bf7a7846 GIT binary patch literal 233 zcmYL@y-EZ@5QOKzK!kgT3BAB=|3VQF10ydu4g-Uxy=mQn-JWH-$8+D|;(G^Q$<%i+ zFmg5!3+kh&O6p=hKiYjip{x3u@TXY~H`RJj5!}ixW)B}OPxU{3_OMpY3Y15Zz}jx9 zgQf6Yj&JNw{3B@_<<*k-m65SuOC7RC8~ri&6uV}trQuK(kXf1XrUfF3czz1^7E>@*!g`kf=O3)CxGb3AOZ#$feZ&AE@lA|DGb33nv8xc8Hzx{2;x_eenx(7s(xZh zYH>z+VopwcO0vF3rIUYra(+sxeriQQYEfoxYF qNqj+RaYnqMNs4}ad}dx|NqoFsLFFwDo80`A(wtN~kX@gFm;nG|RxM8e literal 0 HcmV?d00001 diff --git a/experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-36.pyc b/experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-36.pyc new file mode 100644 index 0000000000000000000000000000000000000000..835a9f72db84c8188e82ceed813fdbd9f3786fde GIT binary patch literal 2774 zcmb7GOOG2j5+>PCJx4R6dHFFGf*hheq*yTmBuL^x5^NG67zYW001X3;wn)wNxYbhF zlr2rVPx7D0;}7_f-;n>18&6s6!I=A=7K`;(QPMn|_z((R{jsV@7THxF^-ir;J|OqQ zahZ^RkxR`2d;m}OB?wM9O~{Z2bjSiWv;u2r2lmhj9EvuUxPkW>;TE@_5N^*|;B$w& zPe@P_W$ry>@^ev|6YyogD4Za(w8Dq z&>8WWu8Ne0VlWVWHIlk6Cux|)!$}gUc$6OYalTbMkNXa21PLf70ppZ2h{h83BE~rl z9Ka?KUDzZQ=p$@_FUF}~~Yd;(9VK=3lJ@G7t2jJ(d*c!RGS zyumk*>NyGF%y3=!Iowk8gf(BHG2Rw$!Wf%xgLJshOGgZl%XgN@;{;;cJv@S?p*0Vd z;MEl%=YT9tA00@H`i!BlK`iyvm1py})od&9>Yu>h^l3A*ru1>m|2=>He5HW=B4__a zg|F@DFitf~$C`~}>_t;;Wz(6qLG-kBgb+lpO+h7+!BZl)VEli6?tYsInRO#Mi&8GK zR@?;b4;XG0U>996cDRa15yt-oC@T2CS@| z@l5++7^krcLk~jDsvwB<3hh`fby=C7ZeAo+XXSY%yBru`vg-hVA}5eBAg!F7*mJr> z&YbQN)4Q;`=T>iF-N8(uhizZhU?%M*aT+B@+D=DusGUg4XeKwnA{!_+QLKY#Q`tm} zC%q8cfV49R;rpN^!yF>kK0R$*+KVw?*$5_`*@z1O!e53%11sdym#o70bPv2U%vEBY z)-El;r+uRZLkphsBK!r`Q6z^V+k+!1;IIHQgDrzCgKdLtgB=A230@BNxZA{WuF-hL z?iuVG>>FG%xMXnI;IhFLgDcDVTb!QJ75l|4 zRX(YlR6*8Y-TIv7`nIF*15pPA5Y}6WQM4P7{>dm(+KN*p zucIDwQJRflIV@@j%SYmFB<{wPo!mu)IW4cDFa?LT)^3?H({?0}9%&}0@XP{x`bZ>! zJwsW>yDiL=7-vD*P>R{0w&WDZQOrrxGw_-WiP0rC61Di z$C(P8$9VMy?kGEs#zH%zEIt)MX*7_r;k&>eagic(2Bin*%JMPrZfJ%Cr3Hpz+BR3V zjdyk6poG{AT?VdsVI)G`HnR63!2oaeqa+E#eK?)$T@Zv->EhR1zOlE&%w;~?0FQ0D z&ETO8+GM+Q_eF+X%n2-r#GFP$5r(=PhQkq`B#5hF`28eG78bb!t7?BNM?n8-w)oOF zKqucraY2a&AumcyS^U}=rmau?TGgpKM%xF4X_zXQB!&z>0HDa66a@bZEO;G?%9*o* zDwop}FQ;?(Q35>%Js>LVqsH9**5PUSi*!6utq)s26sMJY6Se99)nhnIzL)TgMnmmQ z(s%%rFR_s18VK!-CFHSoQ5~r^({!Ne^wRbWQG;dp`X?IJk?w literal 0 HcmV?d00001 diff --git a/experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-37.pyc b/experiments/test_functions/robot_push_14d/__pycache__/push_function.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..c3f12782cc8c0817ec0e91c12c8719f77a0e39ac GIT binary patch literal 2606 zcmZ`)OK%%D5GJ`#JuTbvD{&m7$f3H}h?@44#6}7<2~apidkC-yuw9EQc^zq`a;Zde zc2CML=&eN$_A!_KfnNF_+7nMLdeXh8qG&opt^5MAi=F*AGbD%P@S8`KN|C_x{qMhg zI{@?-Hs+53U?0BZOAwrJ+9f^e(H`?y&+v?%>6tytvnbkF*Y=#pgd5y^OSm~K`T^+Uhp#X4}&CBrIrW- z35`AS zO=WZtf@n4TZ}oY&8;%B{hRS?zPHGqbMYanOT?D=ILO zvb#|nbUVt7`=Y0;K!{*0mcSxvD3(zyf@n}tM~u5P6RUuf)%Ia;(URan5Mzs;)Gn<> zo3AVcPf;(#1_0sb;ZV=;`RD_)Fh1GYr(J$-+8A01oKFeLauL+ z3C&=$05gpZjSY=WjZKX$30ntV7Ph!u$9}ffINI)L>}u?4T+q0naZ%%<#wCqQ`FJ+A z^6|>rUe@-C#ubgL8do)5fL&o(O&WXwB3R_L7f!a6Eob#?C0otbvUO?lMZTCWPK<0r zTIoi*$d~wXT2EJ?r#*498@aa5S*~?JdsDjUV#-pZ&7hByuBGb}JKNM+H!WywCD-1P z1@OI*-b^>sTdA2^X)!IO)pQ}X)AFS9f@NFTZCN}not96dL?(Jz4h zw$|THYv_S3NWqnZ*@eb@5H5%*GN3?WpK4(t!X%ku?oYK#0Z|DB-k>RF@D~k7*{w*l zx}mrQ?arI}&yW8+dAry7@dtd~?VW2>K>!g4ai<#sg=ptb^z++Z^f$yI8qmMNayI#6 z?;MDa(NKeI@Adi-WGwu3k+W0GO>;bQ&`+Uq3;)7t;V*# zRmu!R=aph&1mD;*N3X)JXO2-8v0nIz39v(OuU=-%{77yi(= z=+N{`0X4|8P&U0XTA~2>ElQ+ky~c9p+5O}w7=+5|C(%jh75Z%vX?pWqV81wt(5lbw{F z4duK6%#PvY^e*6x`aR_git14loQkAMaHd$n)s+3pWWR?7dR5qzzDV2pYyR5SCCi%Y80OmnfqDJ+( z_n!Oc`}oe|9=u(xmY)6hPXfDc82@1$dTivc;f;GB%1~z4Xqu8*O$({jwb#hEy3U$$ z$56I%J~x!JXEj}vT;-wUq2!@dP(DgNN(EI^rO%C~uga=|w5Y18hP0&S)I8F%GFOep ziG6(2SZSD<`BsbB>Ej@C4R8D#AdyT6lQETrt~Q3|XU(Rg3K*o6`I|fKhke~uEep;0 zFfHJ{fj9mPM3q0bQX??%nyCqDDVFWT>9`#)H9Dy2m>U+Cw`ZfKAT_>n_e|43zSt@4 zv3EJCu(e_zd%<{B`Z2FOmi9PmHL0meGN-C-mgbd#S=Lte{}WriiMil8cTJl}_Uj;{ z+F&Q%s>ki^oo<-)^*$%R|34&uPx3D$Ux73%&3zHla5UF>l-JUmDL*` z-MY0BJiK>v?QT|`lE7CyBTgrSwXeqK{wpZ z{5VR2?Y@dKFK+25>NU&cw8BBs?)Ne$X>UhO_eQthdX%{>Nxb#t2T^M|tL7bozHV=} zd(G%y!ABfE`8XGZtC`vsQ&Hx!!X`z_coJ$ufB>eMD-u_P@&^?vRi*M+78>@ z#i#MkVx)E&OFD|9P`9?0!Z2(H%a@k!?q2_JDH;wUjivOGcqxhEB-q&LwYb7~N%z=1=-2C_?h+$St55HHCR?KsjV~$QtLAp3?y?+f2 z1qUr4Lmr-~i_Aofs(ag{Dt+{SeqAD(dp~o1osTy$-FbWzaH;@3da6tBm|a>m^oVzrTt%! z`2|S^0wr>e4bGEXAd&An`c38=ByY{Qe&$2QAkI;&Pk{?vuset{e+Qy^GfYBV6DD^j zD@?XWp8&7n>Lrd+XMJI#Jxq3Vq`$?ID@p(h<*Jadi5yEdcWk5i#dsdXFi)8$A!#Mk z;cJa9e7%!8f^iC;oXQhq@^b*NfM0=@sQa;<+B_T70;vkf<|&PBl-<FQ+h9w2!GshnSHn_%I@S=AJXbd}&2?TNC5@nMN<_p7iL-r-- z5#pJ%(e8FNAw?6}dZa7JHvI|Mc z`y_P|E?$;7T{2ad%!WEoKmj{ct`uR3tcqEKWK^)=QT50|SH5)gYwUPHCVJ*NxDA)q z4F)<2f~*wK3E1f}UkQT8J7G7k(XURwd$g}_;)~N z4&R)DLk33=mP5I4{T$_~0#aA`s)*E6^Qxk%Xi-o#`jgrG`m%U}_j*aBH^NqQ`-o4- z3Ebdk)<6c7daO@f^VkJjz&5NKO5ILc?e8(UUCpJVky{ z>}GN0R~TBJwFdm)O&vag2Z9OqSG!R(*jK3A%p3gIyZdy~AencN2tocH^CXQ@voKjs z=02r=n|Z1q;PeDYS#j;&2e*RNwOjYI;)7d14DPMn`e-#P-v4kl&+2O&i94GW(J7E2 z@BuLvzsEW`qJwkdDSi*{#D%Ge3quFSg9B5uM`w?3?#z{$sU<>(#!@!9jtu<+Wk4IA zN<^2+mf@Cx4YOs#td%vWjStaOrYcL*oZS@FM<OnYhwyhpD; zH^-yxUp!MDgUl&+eyo+-`1=1b2|rHy1AhK}ex3_9%ol9WwAB#|@<3#JsP=MgotW65 z*HU;{9V<1zhaC*89`6}e$EG%iao}~y8hy7ODpe10S?Fzcqx$+Xz36&h*Xgd*`x`k2 zTBDxy>&aGBmos0%2J>T?!)5}D%nJvDsHY}yXqI}>(;!zk$54@#@X1Nl2N)#ZZ0M7D zYxL@Kn|;uBO?3avz;e|y>7ybeFmm%sNy3YJmkHaY7m;^`dTj#RQ;lmTJE1+kQcl z-@-e4xP8~~#@xD@ew-J*Nd1`WMMs^Kd$$P<`O-u)dZHN@vIEF zaI=6`HOIc9HuZW}P~)QJK~g}CCp9=()YMVqqoyN$imd5X5?_^4_Y~ZsKM|GDDW(h~ zAdv732vJpF!U{^5NjW8MELFpA?wK<#VvZs;GsMiV0N^w!nJbIce}Fu~4fc9}I8nk^ zuZcead8g`|Fn@eayi+j01?K4N`MBPW>%D$b59_VI=H}JKcyIpg-|+dzcY(C4{7ykd zJzhC8e?zzB@|e4E+{=cOxx8(>mnYXVe>)uB?(0^R6_H!(>&H7$qbd;3PTbHDZm~BH zZpSiER-(%u$j7r{o`)*bSNQQFZaA%I2MSXEK+3ml{}1v#<4ED8kI+le`C`f}5S$Aybi7)M!2%K7FJv6}WR7gA5>v#d>bBIw8$ z4g^6Q7FwO3;EgYV7_d3>v^y50!Gy0cp%trQd8UP*Z~3rB794%%M=u|S^?^O2yzpT< zvquXkF!1^SV`E^lqlF{W znKU`72^pN>UG@e$CLo+h2N!5e$v(pi1j~@I1!(;RV@zL!8WPe)Uj;gnWoPU`W?&mA zIeRAhQCG5GVSz0WjT+1496B^led$jL;*y<}8UEKA2&Fv@A4Slv-^6^fIUFZqQTjoI z3oe}(XFH7JX8FCrpu3x+G^=&nJ%muEO*HpbMaUBb~>c zaCG*_i5%{8;2P2}&6Hga1(;`<+PLs-iov<2(Oahwvqx?VmW<;B#N-ljw;e_ z!d}8k!;+T7k~-rO0~M*iDgI^Z1BO+2Hu(YOW>QlvAGVHP(o5E%>&%+cV#!U9q z=rk)*oS5Zip}*b^vH7g_J4I!Zqy!%&R#EfykP{9ZG zaXFg2FiX7;R@>F!y1m`{|{o7Q8B)Y6t28Ckb_fz^e5Ofz)RfECKz0 zXL3S_X(B3{L)Mg|0itQI5R_*e0a5cqpmy>F&gUkc0r5PfPunWl zYF6?>H$n`3OEx$^6xr`%TLO3G5jF%4z%Ma@xB^1U&4r)jqCZ1kUMBk z=*>+me0WI$)NjgCP32)Z@G@d~m9)AC?f=OI%&Rt@!@-$p4Xc?`5KTEu_s}}BhuniIV%@yj z(Hs|N2TO-@a4l{ey9^`U1pNp7}$T0)FaXviflDWdj>|}eOX>>Ak zH*+w=5HuW;wtM7+LxeqA_%{WA!ZZ=h8UkA4Pr-F$_86wBygTA^Nr>oT{7)D-J_nNH zih*6rticDY*){W`{1zs>nyG2yT2;2j!+_fiZ(pI1@S+JfX%qH*T*awW8N&q^$fD)m zvtZGy%0AB9!&ypSs!8^plM{Wu@7xRf&ZDmj%v7ZBiQ&z}OHU~0G}JSYa(R3*J(-pt z8Tuc9j6zyVi~KtVyb5^#oYH%kPfNn%B^bWKCl&ZNO1kpI@}d6E%WeSi0W@=UJw6=gkAkJX#nPClv_Btu8K*gdTG3r zpM71LQ3)=?iR1ra7>UVp<_ziCIlH7T^Ju{WYBX<~F6nY}e)1nDLbWjT(I9S~IL=wd zKf_Vn){zR(zv*bw$|2Tr(>%tb4IPfnG# zS+Vrrp{`NM+}_UiU{`ia!q83kdVjclv*}B+xW1f+r*5OQ*mP0sqUZ#DGR&d5%-n*^ fe0Q>dNxtt?k^PX6xQIIU z-g6&yAK!UhUa3|~2A;qAw}1SxT{n#XWMTH$$Xv!74M3Ek%)ZeyCAFFsQmb#TlW+B% zb>pU?Y~_4mC}+=Vx+uBIL&-zQL#d#Alzfy5s;Ej|7)@W5RRw8LRaFgXNu5v&NXyDx zGa4uN@l9j3VW#G}HnX#5KB?aTCzJ_cGN!W7)y8KW<$hr{9pzyVKlM91-FrjbS8XdF zfT@f3D&FXG5M_R9B}QQ2H4_umQY_oC({rDAiP1w%&)l@Qyr(v53R2@M_o-zfur47TObO2$m$ivx zzYa2~jdr4~deq(C>9^vc-sc(`&i?-((cTdP?%#p2|1Tt8&#D?Ho0YYzpIpDb8r-{m zZT(hS%}ckg-@bWkJuR(&@X5{V>%omXYpH$b&PS;?+U>Nq!?d8f(Wu|rZTe9d2irpx zre4(6VK``($!WJnad$XKow&OlHr=cJVf%jSwk7d4R_=!Fm9(052!^`b=?le0$+u?E-|kepS?bH z$FC-C&qsc7>Lnh)XC`jqJ+(~Z!JT`?gKPJUxYR2@sXR4I$R^WD2i?(SewJ_VUz3C$1SHxMZ9*P_mYcpIZ-{we0X4wq}>K80*#` z+8pX_T|pT*v6kOX&4&cc1^}gB1XtjqUqbHJrFZ+?C~gn8NAD3s5x?&llS}jSRJ!$O z`CTBk)sNm=-rfz!2u8#HZf7_+s5ltr$WnO3&nzW@l03u37|5AtjF*f}D_cs#%A}gR z(J(btFIe9~qN&cnIy!w_2bb%}=bhTajR-RS6R`F_Cn28y^5 z@L{MSb77>Hz=go>j>6R6i9&s?6}NO%m|S{V$k(n;wjWb3m}Gh%3Op8K?m`1vjF4 zr*>j<<53G3D!`UUFtt&36SLg0eftjeD2Z(|DjEJ(xAzsdYtBn=YL28U7afV+k_;sFTiD^oKDME1z@G)eU3Ru65;xpJzYcrT+Lvk;p8i--O1e|0j zSu`haexn68swqs#olI=kwpy?SGGCaRtFi+bkC;xK&2GQ1Pf8XtabK5_ZTdO7qtuP0 zzbu*)LWh2L)F_B>rXJKxH;kJ8MhAum#!7z&ZS?Pw%oi`GHd`_USY|`vmg9$oDOZZH zF;>N_K@=)j@T7WVp>Tcl%V^(lX`x^s!yrgY0gdoZpZQ7uNa_KnHg$nRVanQgC zJjzW_NKE}|FuwR0o60GtaCy;CXg*>tQ$E#axGL{V*KuD>TZR z{P+6}Yo7{PUk7QFnuUBJsr!fyWa_D5fD;n1q{a2yA6*aD*00}5i+8X8Fu1*b{gbt{ zc<1A_EUVw=NZix3h)#hF0ihQ4ev@@_I0q-h^ZOp&+ySYH148%1gZojlCufdt>)a8U z+aoF)nm5_rIx@8A%7BVH6^nwD&BB&rt8CdSYjqup>-W%9rYcL*oZSq$M<R&b*tT0lb}JbCSzIUa5Q{F!pkWlov9Vy)i5*Z+@6_)$C@@$(<@^IWiDK4*Jo zt&X6N2O`@;t(NKD*u)0Cnivi&NMe2uI~Z64I`Y;NJ2B~dNjzdrzFTi8Rd2x!A9VU* zePf0GYkjEebTI0}&5Q%hP>+Z8cq^>SnXh7lsYg?X%>)*y*BXsN*y;=p&C(!z6l6-} z7%I{dKABg2gh8^+hR#^9Ca*lR*#~XcGwZ|E@?lE1PgYKyO<=%p&^S+}sS?h~CnoM%p5XzOp*+$n&r7i!0GC;0t^>B?=bj68#DU-ik3TBn5an z7CaPFG^Pd5AEKNf;0&C^RThx^Zj5LE8r7Woikj3LTtbbDn!9lUHJ;SqY*AB3jgOk1 z^eM7tP>Fq2O59WMcK%cl(kmv6)Df#c21ZmBsIY<(W>QXw981;kJMq|=78CP{$juxx z3(&-xN-cF|xteE;h~C~9jyuTWdGE4#4Ul*qMb7-uWpO&eJO}3F%$rfY8`THHxZbL_ zVQ;a0burPMU;Hhfe|R5Q`vICFUJVEEIS|xs>aGNT;+<)tzu}}VBZLpKACrifdODkBZMG9ZPrh&<3gVVfWd0CuL?6w7)mfnBu^_I=iBn@>vXiAF)5)8BmBn0{!%iDHk0y#J%=Kc70w>i&0F3F`FEK`7 z8S=IS8?a61`jLel87Kx4A(Og+dCym(5{o{4_coa|RvVhddUOeMHRP@JW{^rtj& z`OfNG7S+c4SbU0{l{RzID1$zLI4>&C#T|n43IdkN64)-~55$PUM35wzyAOj>)2)m)t4Lx9-VZhLKp=y+!O{lcA?!%;51T_sWHUNI zGcU!S8>dZPn*~zm>o^br!ln9KZ$^#tLhBRje0`s?28+1LZxcIm;v!7C&*FWM)E#Lk zlhmOqfMY^5qy^jzYz(1HT$w?$EN!ys99so^+1W(@0$sCc%{mWlP!o-C)|#Ap{<>%h zsQvZtgCO5yp2A3bUU;q@@3P5|WQXKHZ%Iy+7muciH{zWDlvUk#N)=7++>(}2m|N1n zMB$Jn^^A@yQ8i-U{c@rvtG%$o!kt%|Ugmf}U%sfD~E@{YNCUUShKJ26Jr^Mn;)6D(g+=Ly%^P1~fB& z&0>~^`NN$q*)v1atVl68#m&NSqZa_RnPG(+Jh1xsK;P%ys-S3W)_z3a8pLm5GytTEPXV z^%xhPk1c3?di_p*Lb!W`b_QAr)EupHLSXUGa2L*t0$j*tYg@=H`r7<2b#9Okr*?GuTAKMLhC-<=zu=l^hNKplZhMH?$mg|dHlBAA-}~{=Bsz+D|g{7y6`1=|15uOJXpX@vfJ}|g`WSZHT9;2 z7&m0NX*Q9T^AvtyUPijnSbmJz>~l*d~)nwV2CkjYfRc5kkcCDTBD6WKSCWfaYZ!-w8VXar^xIvOjQ{y;&X8emxIwi zVc_TukPKJ6p0&&xT)vuJGtbL!Dd*Kp{pV=c=H(cl7x2yv_SY!Hj4a^^b?`67w2D)y z@+LNh-}D4=i?W_tFxpjRALs4iETu0sB>SF_6MeStiRboRKwls&sYu_G<7=^(oK((f z#4$k1#pzOTDNw^5lcP%U7?2CsW-r17OlSi zTpEDc2{Yp+@u^sCs+7hH+1bA)&8YuY;HUBLD~!ZsC3VJh-ke=hmwB{c0X3T8rAxZf zT*&_}#QiI5c{qxiCy#Sh@Ly{fb#Zh)(>oYF`(F|t z`ghDaq6Rmdnz?Lc@8V2fn-xp{0(Ff_>JE0cN4v6Ha<|-cuMEd4*P6a0%Nr}%eH6Vw wG3TPwMbQcRWSB#95pav;@ZZhuN8h9H%2g|!Q2gxLh1wgnANb$$%l?`F1!}wblmGw# literal 0 HcmV?d00001 diff --git a/experiments/test_functions/robot_push_14d/push_function.py b/experiments/test_functions/robot_push_14d/push_function.py new file mode 100755 index 0000000..6761a79 --- /dev/null +++ b/experiments/test_functions/robot_push_14d/push_function.py @@ -0,0 +1,75 @@ +from experiments.test_functions.robot_push_14d.push_utils import b2WorldInterface, make_base, create_body, end_effector, run_simulation + +import numpy as np + + +class PushReward: + def __init__(self): + + # domain of this function + self.xmin = [-5., -5., -10., -10., 2., 0., -5., -5., -10., -10., 2., 0., -5., -5.] + self.xmax = [5., 5., 10., 10., 30., 2.*np.pi, 5., 5., 10., 10., 30., 2.*np.pi, 5., 5.] + + # starting xy locations for the two objects + self.sxy = (0, 2) + self.sxy2 = (0, -2) + # goal xy locations for the two objects + self.gxy = [4, 3.5] + self.gxy2 = [-4, 3.5] + + @property + def f_max(self): + # maximum value of this function + return np.linalg.norm(np.array(self.gxy) - np.array(self.sxy)) \ + + np.linalg.norm(np.array(self.gxy2) - np.array(self.sxy2)) + @property + def dx(self): + # dimension of the input + return self._dx + + def __call__(self, argv): + # returns the reward of pushing two objects with two robots + rx = float(argv[0]) + ry = float(argv[1]) + xvel = float(argv[2]) + yvel = float(argv[3]) + simu_steps = int(float(argv[4]) * 10) + init_angle = float(argv[5]) + rx2 = float(argv[6]) + ry2 = float(argv[7]) + xvel2 = float(argv[8]) + yvel2 = float(argv[9]) + simu_steps2 = int(float(argv[10]) * 10) + init_angle2 = float(argv[11]) + rtor = float(argv[12]) + rtor2 = float(argv[13]) + + initial_dist = self.f_max + + world = b2WorldInterface(False) + oshape, osize, ofriction, odensity, bfriction, hand_shape, hand_size = \ + 'circle', 1, 0.01, 0.05, 0.01, 'rectangle', (1, 0.3) + + base = make_base(500, 500, world) + body = create_body(base, world, 'rectangle', (0.5, 0.5), ofriction, odensity, self.sxy) + body2 = create_body(base, world, 'circle', 1, ofriction, odensity, self.sxy2) + + robot = end_effector(world, (rx,ry), base, init_angle, hand_shape, hand_size) + robot2 = end_effector(world, (rx2,ry2), base, init_angle2, hand_shape, hand_size) + (ret1, ret2) = run_simulation(world, body, body2, robot, robot2, xvel, yvel, \ + xvel2, yvel2, rtor, rtor2, simu_steps, simu_steps2) + + ret1 = np.linalg.norm(np.array(self.gxy) - ret1) + ret2 = np.linalg.norm(np.array(self.gxy2) - ret2) + return initial_dist - ret1 - ret2 + + +def main(): + f = PushReward() + x = np.random.uniform(f.xmin, f.xmax) + print('Input = {}'.format(x)) + print('Output = {}'.format(f(x))) + + +if __name__ == '__main__': + main() diff --git a/experiments/test_functions/robot_push_14d/push_utils.py b/experiments/test_functions/robot_push_14d/push_utils.py new file mode 100755 index 0000000..d9c3d7e --- /dev/null +++ b/experiments/test_functions/robot_push_14d/push_utils.py @@ -0,0 +1,234 @@ +import numpy as np +import pygame +from Box2D import * +from Box2D.b2 import * + + +class guiWorld: + def __init__(self, fps): + self.SCREEN_WIDTH, self.SCREEN_HEIGHT = 1000, 1000 + self.TARGET_FPS = fps + self.PPM = 10.0 # pixels per meter + self.screen = pygame.display.set_mode((self.SCREEN_WIDTH, self.SCREEN_HEIGHT), 0, 32) + pygame.display.set_caption('push simulator') + self.clock = pygame.time.Clock() + self.screen_origin = b2Vec2(self.SCREEN_WIDTH / (2 * self.PPM), self.SCREEN_HEIGHT / (self.PPM * 2)) + self.colors = { + b2_staticBody: (255, 255, 255, 255), + b2_dynamicBody: (163, 209, 224, 255) + } + + def draw(self, bodies, bg_color=(64, 64, 64, 0)): + def my_draw_polygon(polygon, body, fixture): + vertices = [(self.screen_origin + body.transform * v) * self.PPM for v in polygon.vertices] + vertices = [(v[0], self.SCREEN_HEIGHT - v[1]) for v in vertices] + color = self.colors[body.type] + if body.userData == "obs": + color = (123, 128, 120, 0) + if body.userData == "hand": + color = (174, 136, 218, 0) + + pygame.draw.polygon(self.screen, color, vertices) + + def my_draw_circle(circle, body, fixture): + position = (self.screen_origin + body.transform * circle.pos) * self.PPM + position = (position[0], self.SCREEN_HEIGHT - position[1]) + color = self.colors[body.type] + if body.userData == "hand": + color = (174, 136, 218, 0) + pygame.draw.circle(self.screen, color, [int(x) for x in + position], int(circle.radius * self.PPM)) + + b2PolygonShape.draw = my_draw_polygon + b2CircleShape.draw = my_draw_circle + # draw the world + self.screen.fill(bg_color) + self.clock.tick(self.TARGET_FPS) + for body in bodies: + for fixture in body.fixtures: + fixture.shape.draw(body, fixture) + pygame.display.flip() + + +# this is the interface to pybox2d +class b2WorldInterface: + def __init__(self, do_gui=False): + self.world = b2World(gravity=(0.0, 0.0), doSleep=True) + self.do_gui = do_gui + self.TARGET_FPS = 100 + self.TIME_STEP = 1.0 / self.TARGET_FPS + self.VEL_ITERS, self.POS_ITERS = 10, 10 + self.bodies = [] + + if do_gui: + self.gui_world = guiWorld(self.TARGET_FPS) + # raw_input() + else: + self.gui_world = None + + def initialize_gui(self): + if self.gui_world == None: + self.gui_world = guiWorld(self.TARGET_FPS) + self.do_gui = True + + def stop_gui(self): + self.do_gui = False + + def add_bodies(self, new_bodies): + """ add a single b2Body or list of b2Bodies to the world""" + if type(new_bodies) == list: + self.bodies += new_bodies + else: + self.bodies.append(new_bodies) + + def step(self, show_display=True, idx=0): + self.world.Step(self.TIME_STEP, self.VEL_ITERS, self.POS_ITERS) + if show_display and self.do_gui: + self.gui_world.draw(self.bodies) + + +class end_effector: + def __init__(self, b2world_interface, init_pos, base, init_angle, hand_shape='rectangle', hand_size=(0.3, 1)): + world = b2world_interface.world + self.hand = world.CreateDynamicBody(position=init_pos, angle=init_angle) + self.hand_shape = hand_shape + self.hand_size = hand_size + # forceunit for circle and rect + if hand_shape == 'rectangle': + rshape = b2PolygonShape(box=hand_size) + self.forceunit = 30.0 + elif hand_shape == 'circle': + rshape = b2CircleShape(radius=hand_size) + self.forceunit = 100.0 + elif hand_shape == 'polygon': + rshape = b2PolygonShape(vertices=hand_size) + else: + raise Exception("%s is not a correct shape" % hand_shape) + + self.hand.CreateFixture( + shape=rshape, + density=.1, + friction=.1 + ) + self.hand.userData = "hand" + + friction_joint = world.CreateFrictionJoint( + bodyA=base, + bodyB=self.hand, + maxForce=2, + maxTorque=2, + ) + b2world_interface.add_bodies(self.hand) + + def set_pos(self, pos, angle): + self.hand.position = pos + self.hand.angle = angle + + def apply_wrench(self, rlvel=(0, 0), ravel=0): + + avel = self.hand.angularVelocity + delta_avel = ravel - avel + torque = self.hand.mass * delta_avel * 30.0 + self.hand.ApplyTorque(torque, wake=True) + + lvel = self.hand.linearVelocity + delta_lvel = b2Vec2(rlvel) - b2Vec2(lvel) + force = self.hand.mass * delta_lvel * self.forceunit + self.hand.ApplyForce(force, self.hand.position, wake=True) + + def get_state(self, verbose=False): + state = list(self.hand.position) + [self.hand.angle] + \ + list(self.hand.linearVelocity) + [self.hand.angularVelocity] + if verbose: + print_state = ["%.3f" % x for x in state] + print + "position, velocity: (%s), (%s) " % \ + ((", ").join(print_state[:3]), (", ").join(print_state[3:])) + + return state + + +def create_body(base, b2world_interface, body_shape, body_size, body_friction, body_density, obj_loc): + world = b2world_interface.world + + link = world.CreateDynamicBody(position=obj_loc) + if body_shape == 'rectangle': + linkshape = b2PolygonShape(box=body_size) + elif body_shape == 'circle': + linkshape = b2CircleShape(radius=body_size) + elif body_shape == 'polygon': + linkshape = b2PolygonShape(vertices=body_size) + else: + raise Exception("%s is not a correct shape" % body_shape) + + link.CreateFixture( + shape=linkshape, + density=body_density, + friction=body_friction, + ) + friction_joint = world.CreateFrictionJoint( + bodyA=base, + bodyB=link, + maxForce=5, + maxTorque=2, + ) + + b2world_interface.add_bodies([link]) + return link + + +def make_base(table_width, table_length, b2world_interface): + world = b2world_interface.world + base = world.CreateStaticBody( + position=(0, 0), + shapes=b2PolygonShape(box=(table_length, table_width)), + ) + + b2world_interface.add_bodies([base]) + return base + + +def add_obstacles(b2world_interface, obsverts): + world = b2world_interface.world + obs = [] + for verts in obsverts: + tmp = world.CreateStaticBody( + position=(0, 0), + shapes=b2PolygonShape(vertices=verts), + ) + tmp.userData = "obs" + obs.append(tmp) + + # add boundaries + x, y = sm.wbpolygon.exterior.xy + minx, maxx, miny, maxy = np.min(x), np.max(x), np.min(y), np.max(y) + centers = [(0, miny - 1), (0, maxy + 1), (minx - 1, 0), (maxx + 1, 0)] + boxlen = [(maxx - minx, 0.5), (maxx - minx, 0.5), (0.5, maxy - miny), (0.5, maxy - miny)] + for (pos, blen) in zip(centers, boxlen): + tmp = world.CreateStaticBody( + position=pos, + shapes=b2PolygonShape(box=blen), + ) + obs.append(tmp) + b2world_interface.add_bodies(obs) + + +def run_simulation(world, body, body2, robot, robot2, xvel, yvel, \ + xvel2, yvel2, rtor, rtor2, simulation_steps, + simulation_steps2): + # simulating push with fixed direction pointing from robot location to body location + desired_vel = np.array([xvel, yvel]) + rvel = b2Vec2(desired_vel[0] + np.random.normal(0, 0.01), desired_vel[1] + np.random.normal(0, 0.01)) + + desired_vel2 = np.array([xvel2, yvel2]) + rvel2 = b2Vec2(desired_vel2[0] + np.random.normal(0, 0.01), desired_vel2[1] + np.random.normal(0, 0.01)) + + tmax = np.max([simulation_steps, simulation_steps2]) + for t in range(tmax + 100): + if t < simulation_steps: + robot.apply_wrench(rvel, rtor) + if t < simulation_steps2: + robot2.apply_wrench(rvel2, rtor2) + world.step() + + return (list(body.position), list(body2.position)) diff --git a/experiments/test_functions/robot_push_14d/rover_function.py b/experiments/test_functions/robot_push_14d/rover_function.py new file mode 100755 index 0000000..88ecd19 --- /dev/null +++ b/experiments/test_functions/robot_push_14d/rover_function.py @@ -0,0 +1,253 @@ +from rover_utils import RoverDomain, PointBSpline, ConstObstacleCost, NegGeom, AABoxes, UnionGeom, AdditiveCosts, \ + ConstCost +import numpy as np + + +def create_cost_small(): + c = np.array([[0.11353145, 0.17251116], + [0.4849413, 0.7684513], + [0.38840863, 0.10730809], + [0.32968556, 0.37542275], + [0.64342773, 0.32438415], + [0.42, 0.35], + [0.38745546, 0.0688907], + [0.05771529, 0.1670573], + [0.48750001, 0.67864249], + [0.5294646, 0.66245226], + [0.88495861, 0.76770809], + [0.71132462, 0.46580745], + [0.02038182, 0.32146063], + [0.34077448, 0.70446464], + [0.61490175, 0.79081785], + [0.37367616, 0.6720441], + [0.14711569, 0.57060365], + [0.76084188, 0.65168123], + [0.51038721, 0.78655373], + [0.50396508, 0.90299952], + [0.23763956, 0.38260748], + [0.40169679, 0.72553068], + [0.59670114, 0.08541569], + [0.5514408, 0.62855134], + [0.84606733, 0.94264543], + [0.8, 0.19590218], + [0.39181603, 0.46357532], + [0.44800403, 0.27380116], + [0.5681913, 0.1468706], + [0.37418262, 0.69210095]]) + + l = c - 0.05 + h = c + 0.05 + + r_box = np.array([[0.5, 0.5]]) + r_l = r_box - 0.5 + r_h = r_box + 0.5 + + trees = AABoxes(l, h) + r_box = NegGeom(AABoxes(r_l, r_h)) + obstacles = UnionGeom([trees, r_box]) + + start = np.zeros(2) + 0.05 + goal = np.array([0.95, 0.95]) + + costs = [ConstObstacleCost(obstacles, cost=20.), ConstCost(0.05)] + cost_fn = AdditiveCosts(costs) + return cost_fn, start, goal + + +def create_small_domain(): + cost_fn, start, goal = create_cost_small() + + n_points = 10 + traj = PointBSpline(dim=2, num_points=n_points) + n_params = traj.param_size + domain = RoverDomain(cost_fn, + start=start, + goal=goal, + traj=traj, + s_range=np.array([[-0.1, -0.1], [1.1, 1.1]])) + + return domain + + +def create_cost_large(): + c = np.array([[0.43143755, 0.20876147], + [0.38485367, 0.39183579], + [0.02985961, 0.22328303], + [0.7803707, 0.3447003], + [0.93685657, 0.56297285], + [0.04194252, 0.23598362], + [0.28049582, 0.40984475], + [0.6756053, 0.70939481], + [0.01926493, 0.86972335], + [0.5993437, 0.63347932], + [0.57807619, 0.40180792], + [0.56824287, 0.75486851], + [0.35403502, 0.38591056], + [0.72492026, 0.59969313], + [0.27618746, 0.64322757], + [0.54029566, 0.25492943], + [0.30903526, 0.60166842], + [0.2913432, 0.29636879], + [0.78512072, 0.62340245], + [0.29592116, 0.08400595], + [0.87548394, 0.04877622], + [0.21714791, 0.9607346], + [0.92624074, 0.53441687], + [0.53639253, 0.45127928], + [0.99892031, 0.79537837], + [0.84621631, 0.41891986], + [0.39432819, 0.06768617], + [0.92365693, 0.72217512], + [0.95520914, 0.73956575], + [0.820383, 0.53880139], + [0.22378049, 0.9971974], + [0.34023233, 0.91014706], + [0.64960636, 0.35661133], + [0.29976464, 0.33578931], + [0.43202238, 0.11563227], + [0.66764947, 0.52086962], + [0.45431078, 0.94582745], + [0.12819915, 0.33555344], + [0.19287232, 0.8112075], + [0.61214791, 0.71940626], + [0.4522542, 0.47352186], + [0.95623345, 0.74174186], + [0.17340293, 0.89136853], + [0.04600255, 0.53040724], + [0.42493468, 0.41006649], + [0.37631485, 0.88033853], + [0.66951947, 0.29905739], + [0.4151516, 0.77308712], + [0.55762991, 0.26400156], + [0.6280609, 0.53201974], + [0.92727447, 0.61054975], + [0.93206587, 0.42107549], + [0.63885574, 0.37540613], + [0.15303425, 0.57377797], + [0.8208471, 0.16566631], + [0.14889043, 0.35157346], + [0.71724622, 0.57110725], + [0.32866327, 0.8929578], + [0.74435871, 0.47464421], + [0.9252026, 0.21034329], + [0.57039306, 0.54356078], + [0.56611551, 0.02531317], + [0.84830056, 0.01180542], + [0.51282028, 0.73916524], + [0.58795481, 0.46527371], + [0.83259048, 0.98598188], + [0.00242488, 0.83734691], + [0.72505789, 0.04846931], + [0.07312971, 0.30147979], + [0.55250344, 0.23891255], + [0.51161315, 0.46466442], + [0.802125, 0.93440495], + [0.9157825, 0.32441602], + [0.44927665, 0.53380074], + [0.67708372, 0.67527231], + [0.81868924, 0.88356194], + [0.48228814, 0.88668497], + [0.39805433, 0.99341196], + [0.86671752, 0.79016975], + [0.01115417, 0.6924913], + [0.34272199, 0.89543756], + [0.40721675, 0.86164495], + [0.26317679, 0.37334193], + [0.74446787, 0.84782643], + [0.55560143, 0.46405104], + [0.73567977, 0.12776233], + [0.28080322, 0.26036748], + [0.17507419, 0.95540673], + [0.54233783, 0.1196808], + [0.76670967, 0.88396285], + [0.61297539, 0.79057776], + [0.9344029, 0.86252764], + [0.48746839, 0.74942784], + [0.18657635, 0.58127321], + [0.10377802, 0.71463978], + [0.7771771, 0.01463505], + [0.7635042, 0.45498358], + [0.83345861, 0.34749363], + [0.38273809, 0.51890558], + [0.33887574, 0.82842507], + [0.02073685, 0.41776737], + [0.68754547, 0.96430979], + [0.4704215, 0.92717361], + [0.72666234, 0.63241306], + [0.48494401, 0.72003268], + [0.52601215, 0.81641253], + [0.71426732, 0.47077212], + [0.00258906, 0.30377501], + [0.35495269, 0.98585155], + [0.65507544, 0.03458909], + [0.10550588, 0.62032937], + [0.60259145, 0.87110846], + [0.04959159, 0.535785]]) + + l = c - 0.025 + h = c + 0.025 + + r_box = np.array([[0.5, 0.5]]) + r_l = r_box - 0.5 + r_h = r_box + 0.5 + + trees = AABoxes(l, h) + r_box = NegGeom(AABoxes(r_l, r_h)) + obstacles = UnionGeom([trees, r_box]) + + start = np.zeros(2) + 0.05 + goal = np.array([0.95, 0.95]) + + costs = [ConstObstacleCost(obstacles, cost=20.), ConstCost(0.05)] + cost_fn = AdditiveCosts(costs) + return cost_fn, start, goal + + +def create_large_domain(force_start=False, + force_goal=False, + start_miss_cost=None, + goal_miss_cost=None): + cost_fn, start, goal = create_cost_large() + + n_points = 30 + traj = PointBSpline(dim=2, num_points=n_points) + n_params = traj.param_size + domain = RoverDomain(cost_fn, + start=start, + goal=goal, + traj=traj, + start_miss_cost=start_miss_cost, + goal_miss_cost=goal_miss_cost, + force_start=force_start, + force_goal=force_goal, + s_range=np.array([[-0.1, -0.1], [1.1, 1.1]])) + return domain + + +def main(): + def l2cost(x, point): + return 10 * np.linalg.norm(x - point, 1) + + domain = create_large_domain(force_start=False, + force_goal=False, + start_miss_cost=l2cost, + goal_miss_cost=l2cost) + n_points = domain.traj.npoints + + raw_x_range = np.repeat(domain.s_range, n_points, axis=1) + + from ebo_core.helper import ConstantOffsetFn, NormalizedInputFn + + # maximum value of f + f_max = 5.0 + f = ConstantOffsetFn(domain, f_max) + f = NormalizedInputFn(f, raw_x_range) + x_range = f.get_range() + + x = np.random.uniform(x_range[0], x_range[1]) + print('Input = {}'.format(x)) + print('Output = {}'.format(f(x))) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/experiments/test_functions/robot_push_14d/rover_utils.py b/experiments/test_functions/robot_push_14d/rover_utils.py new file mode 100755 index 0000000..76b2c03 --- /dev/null +++ b/experiments/test_functions/robot_push_14d/rover_utils.py @@ -0,0 +1,338 @@ +from itertools import izip + +import numpy as np +import scipy.interpolate as si + + +class Trajectory: + + def __init__(self): + pass + + def set_params(self, start, goal, params): + raise NotImplemented + + def get_points(self, t): + raise NotImplemented + + @property + def param_size(self): + raise NotImplemented + + +class PointBSpline(Trajectory): + """ + dim : number of dimensions of the state space + num_points : number of internal points used to represent the trajectory. + Note, internal points are not necessarily on the trajectory. + """ + + def __init__(self, dim, num_points): + self.tck = None + self.d = dim + self.npoints = num_points + + """ + Set fit the parameters of the spline from a set of points. If values are given for start or goal, + the start or endpoint of the trajectory will be forces on those points, respectively. + """ + + def set_params(self, params, start=None, goal=None): + + points = params.reshape((-1, self.d)).T + + if start is not None: + points = np.hstack((start[:, None], points)) + + if goal is not None: + points = np.hstack((points, goal[:, None])) + + self.tck, u = si.splprep(points, k=3) + + if start is not None: + for a, sv in izip(self.tck[1], start): + a[0] = sv + + if goal is not None: + for a, gv in izip(self.tck[1], goal): + a[-1] = gv + + def get_points(self, t): + assert self.tck is not None, "Parameters have to be set with set_params() before points can be queried." + return np.vstack(si.splev(t, self.tck)).T + + @property + def param_size(self): + return self.d * self.npoints + + +def simple_rbf(x, point): + return (1 - np.exp(-np.sum(((x - point) / 0.25) ** 2))) + + +class RoverDomain: + """ + Rover domain defined on R^d + cost_fn : vectorized function giving a scalar cost to states + start : a start state for the rover + goal : a goal state + traj : a parameterized trajectory object offering an interface + to interpolate point on the trajectory + s_range : the min and max of the state with s_range[0] in R^d are + the mins and s_range[1] in R^d are the maxs + """ + + def __init__(self, cost_fn, + start, + goal, + traj, + s_range, + start_miss_cost=None, + goal_miss_cost=None, + force_start=True, + force_goal=True, + only_add_start_goal=True, + rnd_stream=None): + self.cost_fn = cost_fn + self.start = start + self.goal = goal + self.traj = traj + self.s_range = s_range + self.rnd_stream = rnd_stream + self.force_start = force_start + self.force_goal = force_goal + + self.goal_miss_cost = goal_miss_cost + self.start_miss_cost = start_miss_cost + + if self.start_miss_cost is None: + self.start_miss_cost = simple_rbf + if self.goal_miss_cost is None: + self.goal_miss_cost = simple_rbf + + if self.rnd_stream is None: + self.rnd_stream = np.random.RandomState(np.random.randint(0, 2 ** 32 - 1)) + + # return the negative cost which need to be optimized + def __call__(self, params, n_samples=1000): + self.set_params(params) + + return -self.estimate_cost(n_samples=n_samples) + + def set_params(self, params): + self.traj.set_params(params + self.rnd_stream.normal(0, 1e-4, params.shape), + self.start if self.force_start else None, + self.goal if self.force_goal else None) + + def estimate_cost(self, n_samples=1000): + # get points on the trajectory + points = self.traj.get_points(np.linspace(0, 1.0, n_samples, endpoint=True)) + # compute cost at each point + costs = self.cost_fn(points) + + # estimate (trapezoidal) the integral of the cost along traj + avg_cost = 0.5 * (costs[:-1] + costs[1:]) + l = np.linalg.norm(points[1:] - points[:-1], axis=1) + total_cost = np.sum(l * avg_cost) + + if not self.force_start: + total_cost += self.start_miss_cost(points[0], self.start) + if not self.force_goal: + total_cost += self.goal_miss_cost(points[-1], self.goal) + + return total_cost + + @property + def input_size(self): + return self.traj.param_size + + +class AABoxes: + def __init__(self, lows, highs): + self.l = lows + self.h = highs + + def contains(self, X): + if X.ndim == 1: + X = X[None, :] + + lX = self.l.T[None, :, :] <= X[:, :, None] + hX = self.h.T[None, :, :] > X[:, :, None] + + return (lX.all(axis=1) & hX.all(axis=1)) + + +class NegGeom: + def __init__(self, geometry): + self.geom = geometry + + def contains(self, X): + return ~self.geom.contains(X) + + +class UnionGeom: + def __init__(self, geometries): + self.geoms = geometries + + def contains(self, X): + return np.any(np.hstack([g.contains(X) for g in self.geoms]), axis=1, keepdims=True) + + +class ConstObstacleCost: + def __init__(self, geometry, cost): + self.geom = geometry + self.c = cost + + def __call__(self, X): + return self.c * self.geom.contains(X) + + +class ConstCost: + def __init__(self, cost): + self.c = cost + + def __call__(self, X): + if X.ndim == 1: + X = X[None, :] + return np.ones((X.shape[0], 1)) * self.c + + +class AdditiveCosts: + def __init__(self, fns): + self.fns = fns + + def __call__(self, X): + return np.sum(np.hstack([f(X) for f in self.fns]), axis=1) + + +class GMCost: + def __init__(self, centers, sigmas, weights=None): + self.c = centers + self.s = sigmas + if self.s.ndim == 1: + self.s = self.s[:, None] + self.w = weights + if weights is None: + self.w = np.ones(centers.shape[0]) + + def __call__(self, X): + if X.ndim == 1: + X = X[None, :] + + return np.exp(-np.sum(((X[:, :, None] - self.c.T[None, :, :]) / self.s.T[None, :, :]) ** 2, axis=1)).dot(self.w) + + +def plot_2d_rover(roverdomain, ngrid_points=100, ntraj_points=100, colormap='RdBu', draw_colorbar=False): + import matplotlib.pyplot as plt + # get a grid of points over the state space + points = [np.linspace(mi, ma, ngrid_points, endpoint=True) for mi, ma in zip(*roverdomain.s_range)] + grid_points = np.meshgrid(*points) + points = np.hstack([g.reshape((-1, 1)) for g in grid_points]) + + # compute the cost at each point on the grid + costs = roverdomain.cost_fn(points) + + # get the cost of the current trajectory + traj_cost = roverdomain.estimate_cost() + + # get points on the current trajectory + traj_points = roverdomain.traj.get_points(np.linspace(0., 1.0, ntraj_points, endpoint=True)) + + # set title to be the total cost + plt.title('traj cost: {0}'.format(traj_cost)) + print('traj cost: {0}'.format(traj_cost)) + # plot cost function + cmesh = plt.pcolormesh(grid_points[0], grid_points[1], costs.reshape((ngrid_points, -1)), cmap=colormap) + if draw_colorbar: + plt.gcf().colorbar(cmesh) + # plot traj + plt.plot(traj_points[:, 0], traj_points[:, 1], 'g') + # plot start and goal + plt.plot([roverdomain.start[0], roverdomain.goal[0]], (roverdomain.start[1], roverdomain.goal[1]), 'ok') + + return cmesh + + +def generate_verts(rectangles): + poly3d = [] + all_faces = [] + vertices = [] + for l, h in zip(rectangles.l, rectangles.h): + verts = [[l[0], l[1], l[2]], [l[0], h[1], l[2]], [h[0], h[1], l[2]], [h[0], l[1], l[2]], + [l[0], l[1], h[2]], [l[0], h[1], h[2]], [h[0], h[1], h[2]], [h[0], l[1], h[2]]] + + faces = [[0, 1, 2, 3], [0, 3, 7, 4], [3, 2, 6, 7], [7, 6, 5, 4], [1, 5, 6, 2], [0, 4, 5, 1]] + + vert_ind = [[0, 1, 2], [0, 2, 3], [0, 3, 4], [4, 3, 7], [7, 3, 2], [2, 6, 7], + [7, 5, 4], [7, 6, 5], [2, 5, 6], [2, 1, 5], [0, 1, 4], [1, 4, 5]] + + plist = [[verts[vert_ind[ix][iy]] for iy in range(len(vert_ind[0]))] for ix in range(len(vert_ind))] + faces = [[verts[faces[ix][iy]] for iy in range(len(faces[0]))] for ix in range(len(faces))] + + poly3d = poly3d + plist + vertices = vertices + verts + all_faces = all_faces + faces + + return poly3d, vertices, all_faces + + +def plot_3d_forest_rover(roverdomain, rectangles, ntraj_points=100): + from matplotlib import pyplot as plt + from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection + + # get the cost of the current trajectory + traj_cost = roverdomain.estimate_cost() + + # get points on the current trajectory + traj_points = roverdomain.traj.get_points(np.linspace(0., 1.0, ntraj_points, endpoint=True)) + + # convert the rectangles into lists of vertices for matplotlib + poly3d, verts, faces = generate_verts(rectangles) + + ax = plt.gcf().add_subplot(111, projection='3d') + + # plot start and goal + ax.scatter((roverdomain.start[0], roverdomain.goal[0]), + (roverdomain.start[1], roverdomain.goal[1]), + (roverdomain.start[2], roverdomain.goal[2]), c='k') + + # plot traj + seg = (zip(traj_points[:-1, :], traj_points[1:, :])) + ax.add_collection3d(Line3DCollection(seg, colors=[(0, 1., 0, 1.)] * len(seg))) + + # plot rectangles + ax.add_collection3d(Poly3DCollection(poly3d, facecolors=(0.7, 0.7, 0.7, 1.), linewidth=0.5)) + + # set limits of axis to be the same as domain + s_range = roverdomain.s_range + ax.set_xlim(s_range[0][0], s_range[1][0]) + ax.set_ylim(s_range[0][1], s_range[1][1]) + ax.set_zlim(s_range[0][2], s_range[1][2]) + + +def main(): + import matplotlib.pyplot as plt + center = np.array([[1., 1.], [1., 0.0]]) + sigma = np.ones(2) * 0.5 + cost_fn = GMCost(center, sigma) + start = np.zeros(2) + 0.1 + goal = np.ones(2) * 1 - 0.1 + + traj = PointBSpline(dim=2, num_points=3) + p = np.array([[0.1, 0.5], [0.3, 1.3], [0.75, 1.2]]) + traj.set_params(start, goal, p.flatten()) + + domain = RoverDomain(cost_fn, + start=start, + goal=goal, + traj=traj, + s_range=np.array([[0., 0.], [2., 2.]])) + + plt.figure() + plot_2d_rover(domain) + plt.plot(p[:, 0], p[:, 1], '*g') + plt.show() + + +if __name__ == "__main__": + main() diff --git a/experiments/test_functions/robot_push_14d/simple_functions.py b/experiments/test_functions/robot_push_14d/simple_functions.py new file mode 100755 index 0000000..bcbf74a --- /dev/null +++ b/experiments/test_functions/robot_push_14d/simple_functions.py @@ -0,0 +1,120 @@ +import os + +import numpy as np + +from ebo_core.bo import global_minimize +from gp_tools.gp import DenseKernelGP +from gp_tools.representation import DenseL1Kernel + + +class SimpleQuadratic(object): + def __init__(self, dx, z, k): + self.dx = dx + self.z = z + self.k = k + + def __call__(self, x): + x = np.squeeze(x) + assert len(x) == self.dx + f = 0 + for i in range(self.dx): + f -= ((x[i] - 0.5) ** 2.0) / self.dx + return f + + +class SampledGpFunc(object): + def __init__(self, x_range, dx, z, k, n, sigma): + self.dx = dx + self.z = z + self.k = k + self.sigma = sigma + self.x_range = x_range + kern = DenseL1Kernel(self.z, self.k) + + X = np.random.uniform(x_range[0], x_range[1], (n, dx)) + kxx = kern(X) + sigma ** 2 * np.eye(X.shape[0]) + y = np.random.multivariate_normal(np.zeros(n), kxx).T + + self.gp = DenseKernelGP(X, y, sigma=sigma, kern=kern) + self.gp.fit() + self.get_max() + + def get_max(self): + x = self.x_range[0].copy() + all_cat = np.unique(self.z) + for a in all_cat: + active = self.z == a + k1 = DenseL1Kernel(self.z[active], self.k[active]) + af = lambda x: np.array(k1(x, self.gp.X[:, active])).dot(self.gp.alpha) + x[active] = np.squeeze(global_minimize(af, self.x_range[:, active], 10000)) + + self.argmax = x + self.f_max = -np.squeeze(np.array(self.gp.kern(x, self.gp.X)).dot(self.gp.alpha)) + + def __call__(self, x): + if x.ndim == 1: + n = 1 + else: + n = x.shape[0] + kXn = np.array(self.gp.kern(x, self.gp.X)) + mu = kXn.dot(self.gp.alpha) + f = mu + f += np.random.normal(size=n) * self.sigma + return -np.squeeze(f) + + +def sample_z(dx): + z = np.zeros(dx, dtype=int) + cnt = 1 + samecnt = 1 + for i in range(1, dx): + if (samecnt < 3 and np.random.rand() < 0.5) or (samecnt < 4 and np.random.rand() < 0.1): + z[i] = z[i - 1] + samecnt += 1 + else: + z[i] = cnt + cnt += 1 + samecnt = 1 + return z + + +def save_sampled_gp_funcs(dx, n=50, nfunc=1, isplot=1, dirnm='mytests'): + import cPickle as pic + for i in range(nfunc): + sigma = 0.01 + z = sample_z(dx) + k = np.array([10] * dx) + x_range = np.matlib.repmat([[0.], [1.]], 1, dx) + f = SampledGpFunc(x_range, dx, z, k, n, sigma) + filenm = os.path.join(dirnm, str(i) + '_' + str(dx) + '_f.pk') + pic.dump(f, open(filenm, 'wb')) + + if isplot: + plot_f(f) + + return f + + +def plot_f(f, filenm='test_function.eps'): + # only for 2D functions + import matplotlib.pyplot as plt + import matplotlib + font = {'size': 20} + matplotlib.rc('font', **font) + + delta = 0.005 + x = np.arange(0.0, 1.0, delta) + y = np.arange(0.0, 1.0, delta) + nx = len(x) + X, Y = np.meshgrid(x, y) + + xx = np.array((X.ravel(), Y.ravel())).T + yy = f(xx) + + plt.figure() + plt.contourf(X, Y, yy.reshape(nx, nx), levels=np.linspace(yy.min(), yy.max(), 40)) + plt.xlim([0, 1]) + plt.ylim([0, 1]) + plt.colorbar() + plt.scatter(f.argmax[0], f.argmax[1], s=180, color='k', marker='+') + plt.savefig(filenm) diff --git a/experiments/test_functions/speed_reducer.py b/experiments/test_functions/speed_reducer.py new file mode 100644 index 0000000..c99a443 --- /dev/null +++ b/experiments/test_functions/speed_reducer.py @@ -0,0 +1,121 @@ +from __future__ import division, print_function +from experiments.test_functions.experiment_configuration import sample_speed_reducer_points +import numpy as np +import torch +class Problem(object): + def __init__(self, dimension, lower_bounds, upper_bounds): + self.dimension = dimension + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + pass + def __call__(self, x_unorder): + x = [] + x.append(2.6 + (x_unorder[1]+1)/2) + x.append(0.7 + ((x_unorder[2]+1) * 0.1)/2) + x.append(x_unorder[0] + 17) + x.append(7.3 + (x_unorder[3]+1)/2) + x.append(7.3 + (x_unorder[4]+1)/2) + x.append(2.9 + (x_unorder[5]+1)/2) + x.append(5 + ((x_unorder[6] + 1) * 0.5)/2) + func_value = 0.7854 * x[0] * (x[1]**2) * (3.333*(x[2]**2) + 14.9334*(x[2]) - 43.0934) \ + - 1.508 * x[0] * (x[5]**2 + x[6]**2) + 7.4777 * (x[5] ** 3 + x[6] ** 3) \ + + 0.7854 * (x[3] * (x[5] ** 2) + x[4] * (x[6] ** 2)) + constraint_cost = 0 + if (27/(x[0] * (x[1] ** 2) * x[2]) >=1): + constraint_cost += 1000 + if (397.5/(x[0] * (x[1] ** 2) * (x[2]**2)) >= 1): + constraint_cost += 1000 + if (1.93 * (x[3] ** 3)/(x[1] * x[2] * (x[5] ** 4)) >= 1): + constraint_cost += 1000 + if (1.93 * (x[4] ** 3)/(x[1] * x[2] * (x[6] ** 4)) >= 1): + constraint_cost += 1000 + if (np.sqrt(((745*x[3]/(x[1]*x[2]))**2) +16.9 * 1e6)/(110*(x[5]**3)) >= 1): + constraint_cost += 1000 + if (np.sqrt(((745*x[4]/(x[1]*x[2]))**2) +157.5 * 1e6)/(85*(x[6]**3)) >= 1): + constraint_cost += 1000 + if (x[1]*x[2]/40 >= 1): + constraint_cost += 1000 + if (5*x[1]/x[0] >=1): + constraint_cost += 1000 + if (x[0]/(12*x[1]) >= 1): + constraint_cost += 1000 + if ((1.5*x[5] + 1.9)/x[3] >= 1): + constraint_cost += 1000 + if ((1.1*x[6] + 1.9)/x[4] >= 1): + constraint_cost += 1000 + return func_value#/100 + + + +class SpeedReducer(object): + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + lower_bounds = [] + upper_bounds = [] + # discrete variables + lower_bounds.append(0) + upper_bounds.append(11) + + + for i in range(6): + lower_bounds.append(-1) + upper_bounds.append(1) + assert len(lower_bounds) == 7 + assert len(upper_bounds) == 7 + print(lower_bounds) + print(upper_bounds) + self.num_discrete = 1 + self.num_continuous = 6 + + self.problem = Problem(dimension=7, lower_bounds=lower_bounds, upper_bounds=upper_bounds) + self.problem.dimension = 7 + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [12] + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.generate_random_points(n_points=10, random_seed=random_seed).float() + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x = x_unorder.numpy() + evaluation = self.problem(x) + print(evaluation) + return torch.tensor(evaluation).float() + + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + random_point.append(torch.randint(0, 12, (1,))) + for i in range(6): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + +if __name__ == '__main__': + mixobj = SpeedReducer(44,'analog_zhi') + print(mixobj.suggested_init) + print(mixobj.suggested_init[0].type()) + print(mixobj.evaluate(mixobj.suggested_init[0])) + print(mixobj.evaluate(mixobj.suggested_init[5])) + diff --git a/experiments/test_functions/weld_design.py b/experiments/test_functions/weld_design.py new file mode 100644 index 0000000..e86a63d --- /dev/null +++ b/experiments/test_functions/weld_design.py @@ -0,0 +1,126 @@ +from __future__ import division, print_function +import numpy as np +import torch +class Problem(object): + def __init__(self, dimension, lower_bounds, upper_bounds): + self.dimension = dimension + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + pass + def __call__(self, x_unorder): + w = x_unorder[0] # w + m = x_unorder[1] # m + h = 0.0625 + ((x_unorder[2] + 1) * (2 - 0.0625))/2 # h + l = ((x_unorder[3] + 1) * (10 - 0))/2# l + t = 2 + ((x_unorder[4] + 1) * (20 - 2))/2# t + b = 0.0625 + ((x_unorder[5] + 1) * (2 - 0.0625))/2 # b + print(f"w: {w}, m: {m}, h: {h}, l: {l}, t: {t}, b: {b}") + if w == 0: + A = np.sqrt(2) * h * l + J = np.sqrt(2) * h * l * (((h+t)**2)/4 + (l**2)/12) + R = np.sqrt(l**2 + (h+t) ** 2)/2 + costheta = l/(2*R) + else: + A = np.sqrt(2) * h * (t + l) + J = np.sqrt(2) * h * l * (((h+t)**2)/4 + (l**2)/12) + np.sqrt(2) * h * t * (((h+l)**2)/4 + (t**2)/12) + R = max(np.sqrt(l**2 + (h+t) ** 2)/2, np.sqrt(t**2 + (h+l) ** 2)/2) + costheta = l/(2*R) + if (m == 0): + C_1, C_2, E, sd, G = 0.1047, 0.0481, 3e7, 3e4 , 12e6 + elif (m == 1): + C_1, C_2, E, sd, G = 0.0489, 0.0224, 14e6, 8e3, 6e6 + elif (m == 2): + C_1, C_2, E, sd, G = 0.5235, 0.2405, 1e7, 5e3, 4e6 + else: + C_1, C_2, E, sd, G = 0.5584, 0.2566, 16e6, 8e3, 6e6 + L = 14 + F = 6000 + delta_max = 0.25 + sigma = 6*F*L/((t**2) * b) + delta = 4*F*(L**3)/(E * (t**3) * b) + pc = 4.013*t*(b**3)*(np.sqrt(E*G))** (1-(t*np.sqrt(E)/4*L*np.sqrt(G)))/(6*(L**2)) + tau1 = F/A + tau2 = F * (L + 0.05*l) * R/J + tau = np.sqrt(tau1**2 + tau2**2 + 2 * tau1*tau2*costheta) + func_value = (1+C_1) * (w*t + l) * (h ** 2) + C_2 * t * b * (14 + l) + constraint_cost = (0.577 * sd < tau) * 10 + (sd < sigma) * 10 + (b < h) * 10 + (pc < F) * 10 + (delta_max < delta) * 10 + constraint_cost = 0 + return func_value + constraint_cost + +class Weld_Design(object): + def __init__(self, random_seed=None, problem_id=None): + self.random_seed = random_seed + self.problem_id = problem_id + lower_bounds = [] + upper_bounds = [] + # discrete variables + lower_bounds.append(0) # weld_config (w) + upper_bounds.append(1) + lower_bounds.append(0) # bulk_material (b) + upper_bounds.append(3) + + for i in range(4): + lower_bounds.append(-1) # weld_thickness (h) + upper_bounds.append(1) + + assert len(lower_bounds) == 6 + assert len(upper_bounds) == 6 + print(lower_bounds) + print(upper_bounds) + self.problem = Problem(dimension=6, lower_bounds=lower_bounds, upper_bounds=upper_bounds) + self.problem.dimension = 6 + self.num_discrete = 2 + self.num_continuous = 4 + print(f"num_discrete: {self.num_discrete}, num_continuous: {self.num_continuous}") + self.n_vertices = [2, 4] + # for i in range(2): + # self.n_vertices.append(10) # m values + self.n_vertices = np.array(self.n_vertices) + self.suggested_init = self.generate_random_points(n_points=10, random_seed=random_seed) + self.adjacency_mat = [] + self.fourier_freq = [] + self.fourier_basis = [] + self.random_seed_info = str(random_seed).zfill(4) + for i in range(len(self.n_vertices)): + n_v = self.n_vertices[i] + #print(n_v) + adjmat = torch.diag(torch.ones(n_v - 1), -1) + torch.diag(torch.ones(n_v - 1), 1) + self.adjacency_mat.append(adjmat) + laplacian = torch.diag(torch.sum(adjmat, dim=0)) - adjmat + eigval, eigvec = torch.symeig(laplacian, eigenvectors=True) + self.fourier_freq.append(eigval) + self.fourier_basis.append(eigvec) + + def evaluate(self, x_unorder): + if x_unorder.dim() == 2: + x_unorder = x_unorder.squeeze(0) + x = x_unorder.numpy().copy() + print(f"evaluating {x}....") + evaluation = self.problem(x) + return torch.tensor(evaluation).float() + def sample_points(self, n_points, random_seed=None): + if random_seed is not None: + rng_state = torch.get_rng_state() + torch.manual_seed(random_seed) + init_points = [] + for _ in range(n_points): + random_point = [] + random_point.append(torch.randint(0, 2, (1,))) + random_point.append(torch.randint(0, 4, (1,))) + + for i in range(4): + random_point.append(torch.FloatTensor(1).uniform_(-1, 1)) + init_points.append(random_point) + return torch.tensor(init_points).float() + + + def generate_random_points(self, n_points, random_seed=None): + return self.sample_points(n_points, random_seed=self.random_seed if random_seed is None else random_seed).float() + + +if __name__ == '__main__': + mixobj = Weld_Design(44,'weld_design') + print(mixobj.suggested_init) + print(mixobj.evaluate(mixobj.suggested_init[0])) + print(mixobj.evaluate(mixobj.suggested_init[5])) + diff --git a/main.py b/main.py new file mode 100644 index 0000000..fc01212 --- /dev/null +++ b/main.py @@ -0,0 +1,188 @@ +import sys +import time +import argparse +import os +import torch + +from GPmodel.kernels.mixeddiffusionkernel import MixedDiffusionKernel +from GPmodel.models.gp_regression import GPRegression +from GPmodel.sampler.sample_mixed_posterior import posterior_sampling + +from GPmodel.sampler.tool_partition import group_input +from GPmodel.inference.inference import Inference + +from acquisition.acquisition_optimization import next_evaluation +from acquisition.acquisition_functions import expected_improvement +from acquisition.acquisition_marginalization import inference_sampling + +from config import experiment_directory +from utils import model_data_filenames, load_model_data, displaying_and_logging + +from experiments.random_seed_config import generate_random_seed_coco +from experiments.test_functions.mixed_integer import MixedIntegerCOCO +from experiments.test_functions.weld_design import Weld_Design +from experiments.test_functions.speed_reducer import SpeedReducer +from experiments.test_functions.pressure_vessel_design import Pressure_Vessel_Design +# from experiments.test_functions.push_robot_14d import Push_robot_14d +from experiments.test_functions.nn_ml_datasets import NN_ML_Datasets +from experiments.test_functions.em_func import EM_func + + + + +def HyBO(objective=None, n_eval=200, path=None, parallel=False, store_data=True, problem_id=None, **kwargs): + """ + :param objective: + :param n_eval: + :param path: + :param parallel: + :param kwargs: + :return: + """ + acquisition_func = expected_improvement + + n_vertices = adj_mat_list = None + eval_inputs = eval_outputs = log_beta = sorted_partition = lengthscales = None + time_list = elapse_list = pred_mean_list = pred_std_list = pred_var_list = None + + if objective is not None: + exp_dir = experiment_directory() + objective_id_list = [objective.__class__.__name__] + if hasattr(objective, 'random_seed_info'): + objective_id_list.append(objective.random_seed_info) + if hasattr(objective, 'data_type'): + objective_id_list.append(objective.data_type) + objective_id_list.append('HyBO') + if problem_id is not None: + objective_id_list.append(problem_id) + objective_name = '_'.join(objective_id_list) + model_filename, data_cfg_filaname, logfile_dir = model_data_filenames(exp_dir=exp_dir, + objective_name=objective_name) + + + n_vertices = objective.n_vertices + adj_mat_list = objective.adjacency_mat + grouped_log_beta = torch.ones(len(objective.fourier_freq)) + log_order_variances = torch.zeros((objective.num_discrete + objective.num_continuous)) + fourier_freq_list = objective.fourier_freq + fourier_basis_list = objective.fourier_basis + suggested_init = objective.suggested_init # suggested_init should be 2d tensor + n_init = suggested_init.size(0) + num_discrete = objective.num_discrete + num_continuous = objective.num_continuous + lengthscales = torch.zeros((num_continuous)) + print("******************* initializing kernel ****************") + kernel = MixedDiffusionKernel(log_order_variances=log_order_variances, grouped_log_beta=grouped_log_beta, fourier_freq_list=fourier_freq_list, + fourier_basis_list=fourier_basis_list, lengthscales=lengthscales, + num_discrete=num_discrete, num_continuous=num_continuous) + surrogate_model = GPRegression(kernel=kernel) + eval_inputs = suggested_init + eval_outputs = torch.zeros(eval_inputs.size(0), 1, device=eval_inputs.device) + for i in range(eval_inputs.size(0)): + eval_outputs[i] = objective.evaluate(eval_inputs[i]) + assert not torch.isnan(eval_outputs).any() + log_beta = eval_outputs.new_zeros(num_discrete) + log_order_variance = torch.zeros((num_discrete + num_continuous)) + sorted_partition = [[m] for m in range(num_discrete)] + lengthscale = torch.zeros((num_continuous)) + + time_list = [time.time()] * n_init + elapse_list = [0] * n_init + pred_mean_list = [0] * n_init + pred_std_list = [0] * n_init + pred_var_list = [0] * n_init + + surrogate_model.init_param(eval_outputs) + print('(%s) Burn-in' % time.strftime('%H:%M:%S', time.localtime())) + sample_posterior = posterior_sampling(surrogate_model, eval_inputs, eval_outputs, n_vertices, adj_mat_list, log_order_variance, + log_beta, lengthscale, sorted_partition, n_sample=1, n_burn=1, n_thin=1) + log_order_variance = sample_posterior[1][0] + log_beta = sample_posterior[2][0] + lengthscale = sample_posterior[3][0] + sorted_partition = sample_posterior[4][0] + print('') + else: + surrogate_model, cfg_data, logfile_dir = load_model_data(path, exp_dir=experiment_directory()) + + for _ in range(n_eval): + start_time = time.time() + reference = torch.min(eval_outputs, dim=0)[0].item() + print('(%s) Sampling' % time.strftime('%H:%M:%S', time.localtime())) + sample_posterior = posterior_sampling(surrogate_model, eval_inputs, eval_outputs, n_vertices, adj_mat_list, log_order_variance, + log_beta, lengthscale, sorted_partition, n_sample=10, n_burn=0, n_thin=1) + hyper_samples, log_order_variance_samples, log_beta_samples, lengthscale_samples, partition_samples, freq_samples, basis_samples, edge_mat_samples = sample_posterior + log_order_variance = log_order_variance_samples[-1] + log_beta = log_beta_samples[-1] + lengthscale = lengthscale_samples[-1] + sorted_partition = partition_samples[-1] + print('\n') + # print(hyper_samples[0]) + # print(log_order_variance) + # print(log_beta) + # print(lengthscale) + # print(sorted_partition) + # print('') + + x_opt = eval_inputs[torch.argmin(eval_outputs)] + inference_samples = inference_sampling(eval_inputs, eval_outputs, n_vertices, + hyper_samples, log_order_variance_samples, log_beta_samples, lengthscale_samples, partition_samples, + freq_samples, basis_samples, num_discrete, num_continuous) + suggestion = next_evaluation(objective, x_opt, eval_inputs, inference_samples, partition_samples, edge_mat_samples, + n_vertices, acquisition_func, reference, parallel) + next_eval, pred_mean, pred_std, pred_var = suggestion + + processing_time = time.time() - start_time + print("next_eval", next_eval) + + eval_inputs = torch.cat([eval_inputs, next_eval.view(1, -1)], 0) + eval_outputs = torch.cat([eval_outputs, objective.evaluate(eval_inputs[-1]).view(1, 1)]) + assert not torch.isnan(eval_outputs).any() + time_list.append(time.time()) + elapse_list.append(processing_time) + pred_mean_list.append(pred_mean.item()) + pred_std_list.append(pred_std.item()) + pred_var_list.append(pred_var.item()) + + displaying_and_logging(logfile_dir, eval_inputs, eval_outputs, pred_mean_list, pred_std_list, pred_var_list, + time_list, elapse_list, hyper_samples, log_beta_samples, lengthscale_samples, log_order_variance_samples, store_data) + print('Optimizing %s with regularization %.2E up to %4d visualization random seed : %s' + % (objective.__class__.__name__, objective.lamda if hasattr(objective, 'lamda') else 0, n_eval, + objective.random_seed_info if hasattr(objective, 'random_seed_info') else 'none')) + + +if __name__ == '__main__': + parser_ = argparse.ArgumentParser( + description='Hybrid Bayesian optimization using additive diffusion kernels') + parser_.add_argument('--n_eval', dest='n_eval', type=int, default=220) + parser_.add_argument('--objective', dest='objective') + parser_.add_argument('--problem_id', dest='problem_id', type=str, default=None) + + args_ = parser_.parse_args() + kwag_ = vars(args_) + objective_ = kwag_['objective'] + print(kwag_) + for i in range(25): + if objective_ == 'coco': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = MixedIntegerCOCO(random_seed_, problem_id=kwag_['problem_id']) + elif objective_ == 'weld_design': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = Weld_Design(random_seed_, problem_id=kwag_['problem_id']) + elif objective_ == 'speed_reducer': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = SpeedReducer(random_seed_, problem_id=kwag_['problem_id']) + elif objective_ == 'pressure_vessel': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = Pressure_Vessel_Design(random_seed_, problem_id=kwag_['problem_id']) + #elif objective_ == 'push_robot': + # random_seed_ = sorted(generate_random_seed_coco())[i] + # kwag_['objective'] = Push_robot_14d(random_seed_, problem_id=kwag_['problem_id']) + elif objective_ == 'em_func': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = EM_func(random_seed_, problem_id=kwag_['problem_id']) + elif objective_ == 'nn_ml_datasets': + random_seed_ = sorted(generate_random_seed_coco())[i] + kwag_['objective'] = NN_ML_Datasets(random_seed_, problem_id=kwag_['problem_id']) + else: + raise NotImplementedError + HyBO(**kwag_) diff --git a/requirement.txt b/requirement.txt new file mode 100644 index 0000000..4d5a63f --- /dev/null +++ b/requirement.txt @@ -0,0 +1,12 @@ +gputil +psutil +simanneal +toposort +torch +torchvision +openml +numpy +pygame +box2d-py +coco: https://github.com/numbbo/coco +cma: https://github.com/CMA-ES/pycma \ No newline at end of file diff --git a/utils.py b/utils.py new file mode 100755 index 0000000..854597b --- /dev/null +++ b/utils.py @@ -0,0 +1,64 @@ +import os +import pickle +import time +from datetime import datetime + +import torch +from config import experiment_directory + + +def model_data_filenames(exp_dir, objective_name): + folder_name = objective_name + '_' + datetime.now().strftime('%Y%m%d-%H:%M:%S:%f') + os.makedirs(os.path.join(exp_dir, folder_name)) + logfile_dir = os.path.join(exp_dir, folder_name, 'log') + os.makedirs(logfile_dir) + model_filename = os.path.join(exp_dir, folder_name, 'model.pt') + cfg_data_filename = os.path.join(exp_dir, folder_name, 'data_config.pkl') + return model_filename, cfg_data_filename, logfile_dir + + +def load_model_data(path, exp_dir=experiment_directory()): + if not os.path.exists(path): + path = os.path.join(exp_dir, path) + logfile_dir = os.path.join(path, 'log') + model_filename = os.path.join(path, 'model.pt') + cfg_data_filename = os.path.join(path, 'data_config.pkl') + + model = torch.load(model_filename) + cfg_data_file = open(cfg_data_filename, 'r') + cfg_data = pickle.load(cfg_data_file) + for key, value in pickle.load(cfg_data_file).iteritems(): + if key != 'logfile_dir': + exec(key + '=value') + cfg_data_file.close() + + return model, cfg_data, logfile_dir + + +def save_model_data(model, model_filename, cfg_data, cfg_data_filename): + torch.save(model, model_filename) + f = open(cfg_data_filename, 'w') + pickle.dump(cfg_data, f) + f.close() + + +def displaying_and_logging(logfile_dir, eval_inputs, eval_outputs, pred_mean_list, pred_std_list, pred_var_list, + time_list, elapse_list, hyper_samples, log_beta_samples, lengthscale_samples, log_order_var_samples, store_data=True): + logfile = open(os.path.join(logfile_dir, str(eval_inputs.size(0)).zfill(4) + '.out'), 'w') + for i in range(eval_inputs.size(0)): + min_val, min_ind = torch.min(eval_outputs[:i + 1], 0) + time_str = time.strftime('%H:%M:%S', time.gmtime(time_list[i])) \ + + '(' + time.strftime('%H:%M:%S', time.gmtime(elapse_list[i])) + ') ' + data_str = ('%3d-th : %+12.4f, mean : %+.4E, std : %.4E, var : %.4E, min : %+8.4f(%3d)' % + (i + 1, eval_outputs.squeeze()[i], + pred_mean_list[i], pred_std_list[i], pred_var_list[i], + min_val.item(), min_ind.item() + 1)) + min_str = ' <==== IMPROVED' if i == min_ind.item() else '' + print(time_str + data_str + min_str) + logfile.writelines(time_str + data_str + min_str + '\n') + logfile.close() + if store_data: + pklfilename = os.path.join(logfile_dir, str(eval_inputs.size(0)).zfill(4) + '.pkl') + torch.save({'inputs': eval_inputs, 'outputs': eval_outputs, 'hyper_samples':hyper_samples, 'log_beta_samples': log_beta_samples, 'lengthscale_samples': lengthscale_samples, 'lov':log_order_var_samples}, pklfilename) + + #torch.save({'inputs': eval_inputs, 'outputs': eval_outputs}, pklfilename)