diff --git a/bloptools/bayesian/__init__.py b/bloptools/bayesian/__init__.py index 822f7be..e49a76f 100644 --- a/bloptools/bayesian/__init__.py +++ b/bloptools/bayesian/__init__.py @@ -1,3 +1,3 @@ from .agent import * # noqa F401 -from .devices import * # noqa F401 -from .objective import * # noqa F401 +from .dofs import * # noqa F401 +from .objectives import * # noqa F401 diff --git a/bloptools/bayesian/acquisition/__init__.py b/bloptools/bayesian/acquisition/__init__.py index 854e63f..a66afbb 100644 --- a/bloptools/bayesian/acquisition/__init__.py +++ b/bloptools/bayesian/acquisition/__init__.py @@ -1,7 +1,6 @@ import os import yaml -from botorch.acquisition.objective import ScalarizedPosteriorTransform from . import analytic, monte_carlo from .analytic import * # noqa F401 @@ -38,7 +37,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb acq_func_name = parse_acq_func_identifier(identifier) acq_func_config = config["upper_confidence_bound"] - if config[acq_func_name]["multitask_only"] and (agent.num_tasks == 1): + if config[acq_func_name]["multitask_only"] and (len(agent.objectives) == 1): raise ValueError(f'Acquisition function "{acq_func_name}" is only for multi-task optimization problems!') # there is probably a better way to structure this @@ -47,7 +46,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, best_f=agent.max_scalarized_objective, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {}} @@ -56,7 +55,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, best_f=agent.max_scalarized_objective, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {}} @@ -65,7 +64,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, best_f=agent.max_scalarized_objective, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {}} @@ -74,7 +73,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, best_f=agent.max_scalarized_objective, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {}} @@ -103,7 +102,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, beta=beta, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {"beta": beta}} @@ -114,7 +113,7 @@ def get_acquisition_function(agent, identifier="qei", return_metadata=True, verb constraint=agent.constraint, model=agent.model, beta=beta, - posterior_transform=ScalarizedPosteriorTransform(weights=agent.objective_weights_torch, offset=0), + posterior_transform=agent.targeting_transform, ) acq_func_meta = {"name": acq_func_name, "args": {"beta": beta}} diff --git a/bloptools/bayesian/acquisition/analytic.py b/bloptools/bayesian/acquisition/analytic.py index c80fe4d..780979e 100644 --- a/bloptools/bayesian/acquisition/analytic.py +++ b/bloptools/bayesian/acquisition/analytic.py @@ -1,58 +1,10 @@ import math -import bluesky.plan_stubs as bps -import bluesky.plans as bp import numpy as np import torch from botorch.acquisition.analytic import LogExpectedImprovement, LogProbabilityOfImprovement, UpperConfidenceBound -def list_scan_with_delay(*args, delay=0, **kwargs): - "Accepts all the normal 'scan' parameters, plus an optional delay." - - def one_nd_step_with_delay(detectors, step, pos_cache): - "This is a copy of bluesky.plan_stubs.one_nd_step with a sleep added." - motors = step.keys() - yield from bps.move_per_step(step, pos_cache) - yield from bps.sleep(delay) - yield from bps.trigger_and_read(list(detectors) + list(motors)) - - kwargs.setdefault("per_step", one_nd_step_with_delay) - uid = yield from bp.list_scan(*args, **kwargs) - return uid - - -def default_acquisition_plan(dofs, inputs, dets, **kwargs): - delay = kwargs.get("delay", 0) - args = [] - for dof, points in zip(dofs, np.atleast_2d(inputs).T): - args.append(dof) - args.append(list(points)) - - uid = yield from list_scan_with_delay(dets, *args, delay=delay) - return uid - - -# def sleepy_acquisition_plan(dofs, inputs, dets): - -# args = [] -# for dof, points in zip(dofs, np.atleast_2d(inputs).T): -# args.append(dof) -# args.append(list(points)) - -# for point in inputs: -# args = [] -# for dof, value in zip(dofs, point): -# args.append(dof) -# args.append(value) - -# yield from bps.mv(*args) -# yield from bps.count([*dets, *dofs]) -# yield from bps.sleep(1) - -# return uid - - class ConstrainedUpperConfidenceBound(UpperConfidenceBound): """Upper confidence bound, but scaled by some constraint. NOTE: Because the UCB can be negative, we constrain it by adjusting the Gaussian quantile. diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index ea00a68..77d1d06 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -3,7 +3,7 @@ import warnings from collections import OrderedDict from collections.abc import Mapping -from typing import Callable, Sequence, Tuple +from typing import Callable, Optional, Sequence, Tuple import bluesky.plan_stubs as bps # noqa F401 import bluesky.plans as bp # noqa F401 @@ -15,6 +15,7 @@ import pandas as pd import scipy as sp import torch +from botorch.acquisition.objective import ScalarizedPosteriorTransform from botorch.models.deterministic import GenericDeterministicModel from botorch.models.model_list_gp_regression import ModelListGP from databroker import Broker @@ -22,10 +23,11 @@ from .. import utils from . import acquisition, models, plotting -from .acquisition import default_acquisition_plan -from .devices import DOF, DOFList from .digestion import default_digestion_function -from .objective import Objective, ObjectiveList +from .dofs import DOF, DOFList +from .objectives import Objective, ObjectiveList +from .plans import default_acquisition_plan +from .transforms import TargetingPosteriorTransform warnings.filterwarnings("ignore", category=botorch.exceptions.warnings.InputDataWarning) @@ -47,6 +49,7 @@ def __init__( tolerate_acquisition_errors=False, sample_center_on_init=False, trigger_delay: float = 0, + train_every: int = 1, ): """ A Bayesian optimization agent. @@ -79,9 +82,9 @@ def __init__( # # below are the behaviors of DOFs of each kind and mode: # - # "read": the agent will read the input on every acquisition (all dofs are always read) - # "move": the agent will try to set and optimize over these (there must be at least one of these) - # "input" means that the agent will use the value to make its posterior + # 'read': the agent will read the input on every acquisition (all dofs are always read) + # 'move': the agent will try to set and optimize over these (there must be at least one of these) + # 'input' means that the agent will use the value to make its posterior # # # not read-only read-only @@ -105,6 +108,7 @@ def __init__( self.tolerate_acquisition_errors = tolerate_acquisition_errors + self.train_every = train_every self.trigger_delay = trigger_delay self.sample_center_on_init = sample_center_on_init @@ -113,7 +117,18 @@ def __init__( self.initialized = False self.a_priori_hypers = None - def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models=True, hypers=None): + self.n_last_trained = 0 + + def tell( + self, + data: Optional[Mapping] = {}, + x: Optional[Mapping] = {}, + y: Optional[Mapping] = {}, + metadata: Optional[Mapping] = {}, + append=True, + train=True, + hypers=None, + ): """ Inform the agent about new inputs and targets for the model. @@ -130,84 +145,172 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= train_models: bool Whether to train the models on construction. hypers: - A dict of hyperparameters for the model to assume a priori. + A dict of hyperparameters for the model to assume a priori, instead of training. """ - new_table = pd.DataFrame({**x, **y, **metadata} if metadata is not None else {**x, **y}) + if not data: + if not x and y: + raise ValueError("Must supply either x and y, or data.") + data = {**x, **y, **metadata} + + data = {k: np.atleast_1d(v) for k, v in data.items()} + unique_field_lengths = {len(v) for v in data.values()} + + if len(unique_field_lengths) > 1: + raise ValueError("All supplies values must be the same length!") + + new_table = pd.DataFrame(data) self.table = pd.concat([self.table, new_table]) if append else new_table self.table.index = np.arange(len(self.table)) - self._update_models(train=train_models, a_priori_hypers=hypers) + for obj in self.objectives: + t0 = ttime.monotonic() - def _update_models(self, train=True, skew_dims=None, a_priori_hypers=None): - skew_dims = skew_dims if skew_dims is not None else self.latent_dim_tuples + cached_hypers = obj.model.state_dict() if hasattr(obj, "model") else None - # if self.initialized: - # cached_hypers = self.hypers + obj.model = self.construct_model(obj) - inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) + if len(obj.model.train_targets) >= 2: + t0 = ttime.monotonic() + self.train_model(obj.model, hypers=(None if train else cached_hypers)) + if self.verbose: + print(f"trained model '{obj.name}' in {1e3*(ttime.monotonic() - t0):.00f} ms") - for i, obj in enumerate(self.objectives): - self.table.loc[:, f"{obj.key}_fitness"] = targets = self._get_objective_targets(i) - train_index = ~np.isnan(targets) + # TODO: should this be per objective? + self.construct_classifier() - if not train_index.sum() >= 2: - raise ValueError("There must be at least two valid data points per objective!") + def train_model(self, model, hypers=None, **kwargs): + """Fit all of the agent's models. All kwargs are passed to `botorch.fit.fit_gpytorch_mll`.""" + if hypers is not None: + model.load_state_dict(hypers) + else: + botorch.fit.fit_gpytorch_mll(gpytorch.mlls.ExactMarginalLogLikelihood(model.likelihood, model), **kwargs) + model.trained = True - train_inputs = torch.tensor(inputs[train_index], dtype=torch.double) - train_targets = torch.tensor(targets[train_index], dtype=torch.double).unsqueeze(-1) # .unsqueeze(0) + def construct_model(self, obj, skew_dims=None): + """ + Construct an untrained model for an objective. + """ - # for constructing the log normal noise prior - # target_snr = 2e2 - # scale = 2e0 - # loc = np.log(1 / target_snr**2) + scale**2 + skew_dims = skew_dims if skew_dims is not None else self.latent_dim_tuples - likelihood = gpytorch.likelihoods.GaussianLikelihood( - noise_constraint=gpytorch.constraints.Interval( - torch.tensor(1e-2).square(), - torch.tensor(1 / obj.min_snr).square(), - ), - # noise_prior=gpytorch.priors.torch_priors.LogNormalPrior(loc=loc, scale=scale), - ) + likelihood = gpytorch.likelihoods.GaussianLikelihood( + noise_constraint=gpytorch.constraints.Interval( + torch.tensor(1e-4).square(), + torch.tensor(1 / obj.min_snr).square(), + ), + # noise_prior=gpytorch.priors.torch_priors.LogNormalPrior(loc=loc, scale=scale), + ) - outcome_transform = botorch.models.transforms.outcome.Standardize(m=1) # , batch_shape=torch.Size((1,))) + outcome_transform = botorch.models.transforms.outcome.Standardize(m=1) # , batch_shape=torch.Size((1,))) - obj.model = models.LatentGP( - train_inputs=train_inputs, - train_targets=train_targets, - likelihood=likelihood, - skew_dims=skew_dims, - input_transform=self._subset_input_transform(active=True), - outcome_transform=outcome_transform, - ) + train_inputs = self.train_inputs(active=True) + train_targets = self.train_targets(obj.name) + + safe = ~(torch.isnan(train_inputs).any(axis=1) | torch.isnan(train_targets).any(axis=1)) + + model = models.LatentGP( + train_inputs=train_inputs[safe], + train_targets=train_targets[safe], + likelihood=likelihood, + skew_dims=skew_dims, + input_transform=self.input_transform, + outcome_transform=outcome_transform, + ) + + model.trained = False + + return model + + def construct_classifier(self, skew_dims=None): + skew_dims = skew_dims if skew_dims is not None else self.latent_dim_tuples dirichlet_likelihood = gpytorch.likelihoods.DirichletClassificationLikelihood( self.all_objectives_valid.long(), learn_additional_noise=True ) self.classifier = models.LatentDirichletClassifier( - train_inputs=torch.tensor(inputs).double(), + train_inputs=self.train_inputs(active=True), train_targets=dirichlet_likelihood.transformed_targets.transpose(-1, -2).double(), skew_dims=skew_dims, likelihood=dirichlet_likelihood, - input_transform=self._subset_input_transform(active=True), + input_transform=self.input_transform, ) - if a_priori_hypers is not None: - self._set_hypers(a_priori_hypers) - else: - self._train_models() - # try: - - # except botorch.exceptions.errors.ModelFittingError: - # if self.initialized: - # self._set_hypers(cached_hypers) - # else: - # raise RuntimeError("Could not fit model on initialization!") - + self.train_model(self.classifier) self.constraint = GenericDeterministicModel(f=lambda x: self.classifier.probabilities(x)[..., -1]) - def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): + # def construct_model(self, obj, skew_dims=None, a_priori_hypers=None): + # ''' + # Construct an untrained model for an objective. + # ''' + # skew_dims = skew_dims if skew_dims is not None else self.latent_dim_tuples + + # inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) + + # for i, obj in enumerate(self.objectives): + # values = self.train_targets(i) + # values = np.where(self.all_objectives_valid, values, np.nan) + + # train_index = ~np.isnan(values) + + # if not train_index.sum() >= 2: + # raise ValueError("There must be at least two valid data points per objective!") + + # train_inputs = torch.tensor(inputs[train_index], dtype=torch.double) + # train_values = torch.tensor(values[train_index], dtype=torch.double).unsqueeze(-1) # .unsqueeze(0) + + # # for constructing the log normal noise prior + # # target_snr = 2e2 + # # scale = 2e0 + # # loc = np.log(1 / target_snr**2) + scale**2 + + # likelihood = gpytorch.likelihoods.GaussianLikelihood( + # noise_constraint=gpytorch.constraints.Interval( + # torch.tensor(1e-4).square(), + # torch.tensor(1 / obj.min_snr).square(), + # ), + # # noise_prior=gpytorch.priors.torch_priors.LogNormalPrior(loc=loc, scale=scale), + # ) + + # outcome_transform = botorch.models.transforms.outcome.Standardize(m=1) # , batch_shape=torch.Size((1,))) + + # obj.model = models.LatentGP( + # train_inputs=train_inputs, + # train_targets=self.t, + # likelihood=likelihood, + # skew_dims=skew_dims, + # input_transform=self.input_transform, + # outcome_transform=outcome_transform, + # ) + + # dirichlet_likelihood = gpytorch.likelihoods.DirichletClassificationLikelihood( + # self.all_objectives_valid.long(), learn_additional_noise=True + # ) + + # self.classifier = models.LatentDirichletClassifier( + # train_inputs=torch.tensor(inputs).double(), + # train_targets=dirichlet_likelihood.transformed_targets.transpose(-1, -2).double(), + # skew_dims=skew_dims, + # likelihood=dirichlet_likelihood, + # input_transform=self._subset_input_transform(active=True), + # ) + + # if a_priori_hypers is not None: + # self._set_hypers(a_priori_hypers) + # else: + # self._train_models() + # # try: + + # # except botorch.exceptions.errors.ModelFittingError: + # # if self.initialized: + # # self._set_hypers(cached_hypers) + # # else: + # # raise RuntimeError('Could not fit model on initialization!') + + # self.constraint = GenericDeterministicModel(f=lambda x: self.classifier.probabilities(x)[..., -1]) + + def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True, upsample=1, **acq_func_kwargs): """Ask the agent for the best point to sample, given an acquisition function. Parameters @@ -229,19 +332,21 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): start_time = ttime.monotonic() if self.verbose: - print(f'finding points with acquisition function "{acq_func_name}" ...') + print(f"finding points with acquisition function '{acq_func_name}' ...") if acq_func_type in ["analytic", "monte_carlo"]: - if not self.initialized: + if not all(hasattr(obj, "model") for obj in self.objectives): raise RuntimeError( - f'Can\'t construct non-trivial acquisition function "{acq_func_identifier}"' + f"Can't construct non-trivial acquisition function '{acq_func_identifier}'" f" (the agent is not initialized!)" ) if acq_func_type == "analytic" and n > 1: raise ValueError("Can't generate multiple design points for analytic acquisition functions.") - acq_func, acq_func_meta = self.get_acquisition_function(identifier=acq_func_identifier, return_metadata=True) + acq_func, acq_func_meta = self.get_acquisition_function( + identifier=acq_func_identifier, return_metadata=True, **acq_func_kwargs + ) NUM_RESTARTS = 8 RAW_SAMPLES = 1024 @@ -255,13 +360,14 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): raw_samples=RAW_SAMPLES, # used for intialization heuristic ) - x = candidates.numpy().astype(float) + # this includes both RO and non-RO DOFs + candidates = candidates.numpy() active_dofs_are_read_only = np.array([dof.read_only for dof in self.dofs.subset(active=True)]) - acq_points = x[..., ~active_dofs_are_read_only] - read_only_X = x[..., active_dofs_are_read_only] - acq_func_meta["read_only_values"] = read_only_X + acq_points = candidates[..., ~active_dofs_are_read_only] + read_only_values = candidates[..., active_dofs_are_read_only] + acq_func_meta["read_only_values"] = read_only_values else: acqf_obj = None @@ -283,20 +389,33 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): raise ValueError() # define dummy acqf objective - acqf_obj = None + acqf_obj = 0 acq_func_meta["duration"] = duration = ttime.monotonic() - start_time if self.verbose: - print( - f"found points {acq_points} with acqf {acq_func_meta['name']} in {duration:.01f} seconds (obj = {acqf_obj})" - ) + print(f"found points {acq_points} in {1e3*duration:.01f} ms (obj = {acqf_obj})") if route and n > 1: routing_index = utils.route(self.dofs.subset(active=True, read_only=False).readback, acq_points) acq_points = acq_points[routing_index] - return acq_points, acq_func_meta + if upsample > 1: + idx = np.arange(len(acq_points)) + upsampled_idx = np.linspace(0, len(idx) - 1, upsample * len(idx) - 1) + acq_points = sp.interpolate.interp1d(idx, acq_points, axis=0)(upsampled_idx) + + res = { + "points": acq_points, + "acq_func": acq_func_meta["name"], + "acq_func_kwargs": acq_func_kwargs, + "duration": acq_func_meta["duration"], + "sequential": sequential, + "upsample": upsample, + "read_only_values": acq_func_meta.get("read_only_values"), + } + + return res def acquire(self, acquisition_inputs): """Acquire and digest according to the self's acquisition and digestion plans. @@ -330,7 +449,7 @@ def acquire(self, acquisition_inputs): # compute the fitness for each objective # for index, entry in products.iterrows(): # for obj in self.objectives: - # products.loc[index, objective["key"]] = getattr(entry, objective["key"]) + # products.loc[index, objective['key']] = getattr(entry, objective['key']) except KeyboardInterrupt as interrupt: raise interrupt @@ -341,30 +460,36 @@ def acquire(self, acquisition_inputs): logging.warning(f"Error in acquisition/digestion: {repr(error)}") products = pd.DataFrame(acquisition_inputs, columns=self.dofs.subset(active=True, read_only=False).names) for obj in self.objectives: - products.loc[:, obj.key] = np.nan + products.loc[:, obj.name] = np.nan if not len(acquisition_inputs) == len(products): raise ValueError("The table returned by the digestion function must be the same length as the sampled inputs!") return products + def load_data(self, data_file, append=True, train=True): + new_table = pd.read_hdf(data_file, key="table") + x = {key: new_table.pop(key).tolist() for key in self.dofs.names} + y = {key: new_table.pop(key).tolist() for key in self.objectives.names} + metadata = new_table.to_dict(orient="list") + self.tell(x=x, y=y, metadata=metadata, append=append, train=train) + def learn( self, - acq_func=None, - n=1, - iterations=1, - upsample=1, - train=True, - data_file=None, + acq_func: str = "qei", + n: int = 1, + iterations: int = 1, + upsample: int = 1, + train: bool = True, + append: bool = True, hypers_file=None, - append=True, ): """This returns a Bluesky plan which iterates the learning algorithm, looping over ask -> acquire -> tell. For example: - RE(agent.learn("qr", n=16)) - RE(agent.learn("qei", n=4, iterations=4)) + RE(agent.learn('qr', n=16)) + RE(agent.learn('qei', n=4, iterations=4)) Parameters ---------- @@ -385,32 +510,22 @@ def learn( hyperparameters a priori for the rest of the run, and not try to fit a model. """ - if data_file is not None: - new_table = pd.read_hdf(data_file, key="table") + if self.sample_center_on_init and not self.initialized: + center_inputs = np.atleast_2d(self.dofs.subset(active=True, read_only=False).limits.mean(axis=1)) + new_table = yield from self.acquire(center_inputs) + new_table.loc[:, "acq_func"] = "sample_center_on_init" - elif acq_func is not None: - if self.sample_center_on_init and not self.initialized: - center_inputs = np.atleast_2d(self.dofs.subset(active=True, read_only=False).limits.mean(axis=1)) - new_table = yield from self.acquire(center_inputs) - new_table.loc[:, "acq_func"] = "sample_center_on_init" + for i in range(iterations): + print(f"running iteration {i + 1} / {iterations}") + for single_acq_func in np.atleast_1d(acq_func): + res = self.ask(n=n, acq_func_identifier=single_acq_func, upsample=upsample) + new_table = yield from self.acquire(res["points"]) + new_table.loc[:, "acq_func"] = res["acq_func"] - for i in range(iterations): - print(f"running iteration {i + 1} / {iterations}") - for single_acq_func in np.atleast_1d(acq_func): - x, acq_func_meta = self.ask(n=n, acq_func_identifier=single_acq_func) - new_table = yield from self.acquire(x) - new_table.loc[:, "acq_func"] = acq_func_meta["name"] - - else: - raise ValueError("You must supply either an acquisition function or a filepath!") - - x = {key: new_table.pop(key).tolist() for key in self.dofs.names} - y = {key: new_table.pop(key).tolist() for key in self.objectives.keys} - metadata = new_table.to_dict(orient="list") - - self.tell(x=x, y=y, metadata=metadata, append=append, train_models=True) - - self.initialized = True + x = {key: new_table.pop(key).tolist() for key in self.dofs.names} + y = {key: new_table.pop(key).tolist() for key in self.objectives.names} + metadata = new_table.to_dict(orient="list") + self.tell(x=x, y=y, metadata=metadata, append=append, train=train) def get_acquisition_function(self, identifier, return_metadata=False): """Returns a BoTorch acquisition function for a given identifier. Acquisition functions can be @@ -421,7 +536,12 @@ def get_acquisition_function(self, identifier, return_metadata=False): def reset(self): """Reset the agent.""" self.table = pd.DataFrame() - self.initialized = False + + for obj in self.objectives: + if hasattr(obj, "model"): + del obj.model + + self.n_last_trained = 0 def benchmark( self, output_dir="./", runs=16, n_init=64, learning_kwargs_list=[{"acq_func": "qei", "n": 4, "iterations": 16}] @@ -431,25 +551,16 @@ def benchmark( Parameters ---------- output_dir : - Where to save the optimized agents + Where to save the agent output. runs : int How many benchmarks to run - n_init : int - How many points to sample on reseting the agent. learning_kwargs_list: - A list of kwargs which the agent will run sequentially for each run. + A list of kwargs to pass to the learn method which the agent will run sequentially for each run. """ - # cache_limits = {dof.name: dof.limits for dof in self.dofs} for run in range(runs): - # for dof in self.dofs: - # offset = 0.25 * np.ptp(dof.limits) * np.random.uniform(low=-1, high=1) - # dof.limits = (cache_limits[dof.name][0] + offset, cache_limits[dof.name][1] + offset) - self.reset() - yield from self.learn("qr", n=n_init) - for kwargs in learning_kwargs_list: yield from self.learn(**kwargs) @@ -460,43 +571,40 @@ def model(self): """A model encompassing all the objectives. A single GP in the single-objective case, or a model list.""" return ModelListGP(*[obj.model for obj in self.objectives]) if len(self.objectives) > 1 else self.objectives[0].model + def posterior(self, x): + """A model encompassing all the objectives. A single GP in the single-objective case, or a model list.""" + return self.model.posterior(x) + @property def objective_weights_torch(self): return torch.tensor(self.objectives.weights, dtype=torch.double) - def _get_objective_targets(self, i): - """Returns the targets (what we fit to) for an objective, given the objective index.""" - obj = self.objectives[i] - - targets = self.table.loc[:, obj.key].values.copy() - - # check that targets values are inside acceptable values - valid = (targets > obj.limits[0]) & (targets < obj.limits[1]) - targets = np.where(valid, targets, np.nan) - - # transform if needed - if obj.log: - targets = np.where(valid, np.log(targets), np.nan) - - if obj.minimize: - targets *= -1 - - return targets + @property + def scalarizing_transform(self): + return ScalarizedPosteriorTransform(weights=self.objective_weights_torch, offset=0) @property - def n_objs(self): - """Returns a (num_objectives x n_observations) array of objectives""" - return len(self.objectives) + def targeting_transform(self): + return TargetingPosteriorTransform(weights=self.objective_weights_torch, targets=self.objectives.targets) @property - def objectives_targets(self): - """Returns a (num_objectives x n_obs) array of objectives""" - return torch.cat([torch.tensor(self._get_objective_targets(i))[..., None] for i in range(self.n_objs)], dim=1) + def pseudo_targets(self): + """Targets for the posterior transform""" + return torch.tensor( + [ + self.train_targets(active=True)[..., i].max() + if t == "max" + else self.train_targets(active=True)[..., i].min() + if t == "min" + else t + for i, t in enumerate(self.objectives.targets) + ] + ) @property def scalarized_objectives(self): """Returns a (n_obs,) array of scalarized objectives""" - return (self.objectives_targets * self.objectives.weights).sum(axis=-1) + return self.targeting_transform.evaluate(self.train_targets(active=True)).sum(axis=-1) @property def max_scalarized_objective(self): @@ -561,6 +669,13 @@ def latent_dim_tuples(self): return [tuple(np.where(uinv == i)[0]) for i in range(len(u))] + @property + def input_transform(self): + """ + A bounding transform for all the active DOFs. This is used for model fitting. + """ + return self._subset_input_transform(active=True) + def _subset_input_transform(self, active=None, read_only=None, tags=[]): # torch likes limits to be (2, n_dof) and not (n_dof, 2) torch_limits = torch.tensor(self.dofs.subset(active, read_only, tags).limits.T, dtype=torch.double) @@ -590,7 +705,7 @@ def forget(self, index, train=True): Make the agent forget some index of the data table. """ self.table.drop(index=index, inplace=True) - self._update_models(train=train) + self._construct_models(train=train) def forget_last_n(self, n, train=True): """ @@ -610,7 +725,7 @@ def sampler(self, n, d): def _set_hypers(self, hypers): for obj in self.objectives: - obj.model.load_state_dict(hypers[obj.key]) + obj.model.load_state_dict(hypers[obj.name]) self.classifier.load_state_dict(hypers["classifier"]) @property @@ -620,9 +735,9 @@ def hypers(self): for key, value in self.classifier.state_dict().items(): hypers["classifier"][key] = value for obj in self.objectives: - hypers[obj.key] = {} + hypers[obj.name] = {} for key, value in obj.model.state_dict().items(): - hypers[obj.key][key] = value + hypers[obj.name][key] = value return hypers @@ -658,23 +773,63 @@ def _train_models(self, **kwargs): if self.verbose: print(f"trained models in {ttime.monotonic() - t0:.01f} seconds") + self.n_last_trained = len(self.table) + @property def all_acq_funcs(self): """Description and identifiers for all supported acquisition functions.""" entries = [] for k, d in acquisition.config.items(): ret = "" - ret += f'{d["pretty_name"].upper()} (identifiers: {d["identifiers"]})\n' - ret += f'-> {d["description"]}' + ret += f"{d['pretty_name'].upper()} (identifiers: {d['identifiers']})\n" + ret += f"-> {d['description']}" entries.append(ret) print("\n\n".join(entries)) @property def inputs(self): - """A two-dimensional array of all DOF values.""" + """A DataFrame of all DOF values.""" return self.table.loc[:, self.dofs.names].astype(float) + def train_inputs(self, dof_name=None, **subset_kwargs): + """A two-dimensional tensor of all DOF values.""" + + if dof_name is None: + return torch.cat([self.train_inputs(dof.name) for dof in self.dofs.subset(**subset_kwargs)], dim=-1) + + dof = self.dofs[dof_name] + inputs = self.table.loc[:, dof.name].values.copy() + + # check that inputs values are inside acceptable values + valid = (inputs >= dof.limits[0]) & (inputs <= dof.limits[1]) + inputs = np.where(valid, inputs, np.nan) + + # transform if needed + if dof.log: + inputs = np.where(inputs > 0, np.log(inputs), np.nan) + + return torch.tensor(inputs, dtype=torch.double).unsqueeze(-1) + + def train_targets(self, obj_name=None, **subset_kwargs): + """Returns the values associated with an objective name.""" + + if obj_name is None: + return torch.cat([self.train_targets(obj.name) for obj in self.objectives], dim=-1) + + obj = self.objectives[obj_name] + targets = self.table.loc[:, obj.name].values.copy() + + # check that targets values are inside acceptable values + valid = (targets >= obj.limits[0]) & (targets <= obj.limits[1]) + targets = np.where(valid, targets, np.nan) + + # transform if needed + if obj.log: + targets = np.where(targets > 0, np.log(targets), np.nan) + + return torch.tensor(targets, dtype=torch.double).unsqueeze(-1) + @property def active_inputs(self): """A two-dimensional array of all inputs for model fitting.""" diff --git a/bloptools/bayesian/devices.py b/bloptools/bayesian/dofs.py similarity index 73% rename from bloptools/bayesian/devices.py rename to bloptools/bayesian/dofs.py index d39bc1f..ad3f39f 100644 --- a/bloptools/bayesian/devices.py +++ b/bloptools/bayesian/dofs.py @@ -1,6 +1,7 @@ import time as ttime import uuid from collections.abc import Sequence +from dataclasses import dataclass, field from typing import Tuple, Union import numpy as np @@ -18,17 +19,18 @@ class ReadOnlyError(Exception): def _validate_dofs(dofs): - """Check that a list of DOFs can be combined into a DOFList.""" + dof_names = [dof.name for dof in dofs] # check that dof names are unique - unique_dof_names, counts = np.unique([dof.name for dof in dofs], return_counts=True) + unique_dof_names, counts = np.unique(dof_names, return_counts=True) duplicate_dof_names = unique_dof_names[counts > 1] if len(duplicate_dof_names) > 0: - raise ValueError(f'Duplicate name(s) in supplied dofs: "{duplicate_dof_names}"') + raise ValueError(f"Duplicate name(s) in supplied dofs: {duplicate_dof_names}") return list(dofs) +@dataclass class DOF: """A degree of freedom (DOF), to be used by an agent. @@ -59,41 +61,44 @@ class DOF: DOF will be modeled independently. """ - def __init__( - self, - device: Signal = None, - name: str = None, - description: str = "", - limits: Tuple[float, float] = (-10.0, 10.0), - units: str = None, - read_only: bool = False, - active: bool = True, - tags: list = [], - latent_group=None, - ): + device: Signal = None + description: str = None + name: str = None + limits: Tuple[float, float] = (-10.0, 10.0) + units: str = "" + read_only: bool = False + active: bool = True + tags: list = field(default_factory=list) + log: bool = False + latent_group: str = None + + # Some post-processing. This is specific to dataclasses + def __post_init__(self): self.uuid = str(uuid.uuid4()) - self.description = description - self.name = name if name is not None else device.name if hasattr(device, "name") else self.uuid + if self.name is None: + self.name = self.device.name if hasattr(self.device, "name") else self.uuid - self.device = device if device is not None else Signal(name=self.name) - self.limits = limits - self.read_only = read_only if read_only is not None else True if isinstance(device, SignalRO) else False - self.units = units - self.tags = tags - self.active = active - self.latent_group = latent_group if latent_group is not None else str(uuid.uuid4()) + if self.device is None: + self.device = Signal(name=self.name) + if not self.read_only: + # check that the device has a put method + if isinstance(self.device, SignalRO): + raise ValueError("Must specify read_only=True for a read-only device!") + + if self.latent_group is None: + self.latent_group = str(uuid.uuid4()) + + # all dof degrees of freedom are hinted self.device.kind = "hinted" @property def lower_limit(self): - """The lower limit of the DOF.""" return float(self.limits[0]) @property def upper_limit(self): - """The upper limit of the DOF.""" return float(self.limits[1]) @property @@ -102,7 +107,6 @@ def readback(self): @property def summary(self) -> pd.Series: - """A pandas Series representing the current state of the DOF.""" series = pd.Series(index=DOF_FIELDS) for attr in series.index: series[attr] = getattr(self, attr) @@ -110,32 +114,35 @@ def summary(self) -> pd.Series: @property def label(self) -> str: - """A formal label for plotting.""" - return f"{self.name}{f' [{self.units}]' if self.units is not None else ''}" + return f"{self.name}{f' [{self.units}]' if len(self.units) > 0 else ''}" + @property + def has_model(self): + return hasattr(self, "model") -class DOFList(Sequence): - """A class for handling a list of DOFs.""" +class DOFList(Sequence): def __init__(self, dofs: list = []): _validate_dofs(dofs) self.dofs = dofs - def __getitem__(self, index): - """Get a DOF either by name or its position in the list.""" - if type(index) is str: - return self.dofs[self.names.index(index)] - if type(index) is int: - return self.dofs[index] + def __getitem__(self, i): + if type(i) is int: + return self.dofs[i] + elif type(i) is str: + return self.dofs[self.names.index(i)] + else: + raise ValueError(f"Invalid index {i}. A DOFList must be indexed by either an integer or a string.") def __len__(self): - """Number of DOFs in the list.""" return len(self.dofs) def __repr__(self): - """A table showing the state of each DOF.""" return self.summary.__repr__() + # def _repr_html_(self): + # return self.summary._repr_html_() + @property def summary(self) -> pd.DataFrame: table = pd.DataFrame(columns=DOF_FIELDS) @@ -160,6 +167,10 @@ def names(self) -> list: def devices(self) -> list: return [dof.device for dof in self.dofs] + @property + def device_names(self) -> list: + return [dof.device.name for dof in self.dofs] + @property def lower_limits(self) -> np.array: return np.array([dof.lower_limit for dof in self.dofs]) @@ -202,26 +213,17 @@ def subset(self, active=None, read_only=None, tags=[]): return DOFList([dof for dof, m in zip(self.dofs, self._dof_mask(active, read_only, tags)) if m]) def activate(self, read_only=None, active=None, tags=[]): - """Activate all degrees of freedom with a given tag, active status or read-only status. - - For example, `dofs.activate(tag='kb')` will turn off all dofs which contain the tag 'kb'. - """ for dof in self._subset_dofs(read_only, active, tags): dof.active = True def deactivate(self, read_only=None, active=None, tags=[]): - """The same as .activate(), only in reverse.""" for dof in self._subset_dofs(read_only, active, tags): dof.active = False class BrownianMotion(SignalRO): - """Read-only degree of freedom simulating Brownian motion. - - Parameters - ---------- - theta : float - Determines the autocorrelation of the process; smaller values correspond to faster variation. + """ + Read-only degree of freedom simulating brownian motion """ def __init__(self, name=None, theta=0.95, *args, **kwargs): diff --git a/bloptools/bayesian/kernels.py b/bloptools/bayesian/kernels.py index dc46438..dc9f642 100644 --- a/bloptools/bayesian/kernels.py +++ b/bloptools/bayesian/kernels.py @@ -12,7 +12,7 @@ def __init__( self, num_inputs=1, skew_dims=True, - diag_prior=True, + priors=True, scale_output=True, **kwargs, ): @@ -63,21 +63,23 @@ def __init__( self.n_skew_entries = len(self.skew_matrix_indices[0]) - diag_entries_constraint = gpytorch.constraints.Positive() - raw_diag_entries_initial = ( - diag_entries_constraint.inverse_transform(torch.tensor(1e-1)) + lengthscale_constraint = gpytorch.constraints.Positive() + raw_lengthscale_entries_initial = ( + lengthscale_constraint.inverse_transform(torch.tensor(1e-1)) * torch.ones(self.num_outputs, self.num_inputs).double() ) - self.register_parameter(name="raw_diag_entries", parameter=torch.nn.Parameter(raw_diag_entries_initial)) - self.register_constraint(param_name="raw_diag_entries", constraint=diag_entries_constraint) + self.register_parameter( + name="raw_lengthscale_entries", parameter=torch.nn.Parameter(raw_lengthscale_entries_initial) + ) + self.register_constraint(param_name="raw_lengthscale_entries", constraint=lengthscale_constraint) - if diag_prior: + if priors: self.register_prior( - name="diag_entries_prior", + name="lengthscale_entries_prior", prior=gpytorch.priors.GammaPrior(concentration=2, rate=1), - param_or_closure=lambda m: m.diag_entries, - setting_closure=lambda m, v: m._set_diag_entries(v), + param_or_closure=lambda m: m.lengthscale_entries, + setting_closure=lambda m, v: m._set_lengthscale_entries(v), ) if self.n_skew_entries > 0: @@ -105,8 +107,8 @@ def __init__( ) @property - def diag_entries(self): - return self.raw_diag_entries_constraint.transform(self.raw_diag_entries) + def lengthscale_entries(self): + return self.raw_lengthscale_entries_constraint.transform(self.raw_lengthscale_entries) @property def skew_entries(self): @@ -116,9 +118,9 @@ def skew_entries(self): def outputscale(self): return self.raw_outputscale_constraint.transform(self.raw_outputscale) - @diag_entries.setter - def diag_entries(self, value): - self._set_diag_entries(value) + @lengthscale_entries.setter + def lengthscale_entries(self, value): + self._set_lengthscale_entries(value) @skew_entries.setter def skew_entries(self, value): @@ -128,10 +130,10 @@ def skew_entries(self, value): def outputscale(self, value): self._set_outputscale(value) - def _set_diag_entries(self, value): + def _set_lengthscale_entries(self, value): if not torch.is_tensor(value): - value = torch.as_tensor(value).to(self.raw_diag_entries) - self.initialize(raw_diag_entries=self.raw_diag_entries_constraint.inverse_transform(value)) + value = torch.as_tensor(value).to(self.raw_lengthscale_entries) + self.initialize(raw_lengthscale_entries=self.raw_lengthscale_entries_constraint.inverse_transform(value)) def _set_skew_entries(self, value): if not torch.is_tensor(value): @@ -155,7 +157,7 @@ def skew_matrix(self): @property def diag_matrix(self): D = torch.zeros((self.num_outputs, self.num_inputs, self.num_inputs), dtype=torch.float64) - D[self.diag_matrix_indices] = self.diag_entries.ravel() + D[self.diag_matrix_indices] = self.lengthscale_entries.ravel() ** -1 return D @property diff --git a/bloptools/bayesian/models.py b/bloptools/bayesian/models.py index eff9a64..a50e84e 100644 --- a/bloptools/bayesian/models.py +++ b/bloptools/bayesian/models.py @@ -15,7 +15,23 @@ def __init__(self, train_inputs, train_targets, skew_dims=True, *args, **kwargs) num_inputs=train_inputs.shape[-1], num_outputs=train_targets.shape[-1], skew_dims=skew_dims, - diag_prior=True, + priors=True, + scale=True, + **kwargs + ) + + +class TargetingGP(botorch.models.gp_regression.SingleTaskGP): + def __init__(self, train_inputs, train_targets, skew_dims=True, *args, **kwargs): + super().__init__(train_inputs, train_targets, *args, **kwargs) + + self.mean_module = gpytorch.means.ConstantMean(constant_prior=gpytorch.priors.NormalPrior(loc=0, scale=1)) + + self.covar_module = kernels.LatentKernel( + num_inputs=train_inputs.shape[-1], + num_outputs=train_targets.shape[-1], + skew_dims=skew_dims, + priors=True, scale=True, **kwargs ) diff --git a/bloptools/bayesian/objective.py b/bloptools/bayesian/objective.py deleted file mode 100644 index 47b963f..0000000 --- a/bloptools/bayesian/objective.py +++ /dev/null @@ -1,145 +0,0 @@ -from collections.abc import Sequence -from typing import Tuple, Union - -import numpy as np -import pandas as pd - -numeric = Union[float, int] - -DEFAULT_MINIMUM_SNR = 1e1 -OBJ_FIELDS = ["description", "key", "limits", "weight", "minimize", "log"] - - -class DuplicateKeyError(ValueError): - ... - - -def _validate_objectives(objectives): - keys = [obj.key for obj in objectives] - unique_keys, counts = np.unique(keys, return_counts=True) - duplicate_keys = unique_keys[counts > 1] - if len(duplicate_keys) > 0: - raise DuplicateKeyError(f'Duplicate key(s) in supplied objectives: "{duplicate_keys}"') - - -class Objective: - """A degree of freedom (DOF), to be used by an agent. - - Parameters - ---------- - description: str - The description of the DOF. This is used as a key. - description: str - A longer description for the DOF. - device: Signal, optional - An ophyd device. If None, a dummy ophyd device is generated. - limits: tuple, optional - A tuple of the lower and upper limit of the DOF. If the DOF is not read-only, the agent - will not explore outside the limits. If the DOF is read-only, the agent will reject all - sampled data where the DOF is outside the limits. - read_only: bool - If True, the agent will not try to set the DOF. Must be set to True if the supplied ophyd - device is read-only. - active: bool - If True, the agent will try to use the DOF in its optimization. If False, the agent will - still read the DOF but not include it any model or acquisition function. - units: str - The units of the DOF (e.g. mm or deg). This is only for plotting and general housekeeping. - tags: list - A list of tags. These make it easier to subset large groups of dofs. - latent_group: optional - An agent will fit latent dimensions to all DOFs with the same latent_group. If None, the - DOF will be modeled independently. - """ - - def __init__( - self, - key: str, - description: str = "", - minimize: bool = False, - log: bool = False, - weight: numeric = 1.0, - limits: Tuple[numeric, numeric] = None, - min_snr: numeric = DEFAULT_MINIMUM_SNR, - ): - self.description = description if description is not None else key - self.key = key - self.minimize = minimize - self.log = log - self.weight = weight - self.min_snr = min_snr - - if limits is not None: - self.limits = limits - elif self.log: - self.limits = (0, np.inf) - else: - self.limits = (-np.inf, np.inf) - - @property - def label(self): - return f"{'neg ' if self.minimize else ''}{'log ' if self.log else ''}{self.description}" - - @property - def summary(self): - series = pd.Series() - for col in OBJ_FIELDS: - series[col] = getattr(self, col) - return series - - def __repr__(self): - return self.summary.__repr__() - - @property - def noise(self): - return self.model.likelihood.noise.item() if hasattr(self, "model") else None - - -class ObjectiveList(Sequence): - def __init__(self, objectives: list = []): - _validate_objectives(objectives) - self.objectives = objectives - - def __getitem__(self, i): - return self.objectives[i] - - def __len__(self): - return len(self.objectives) - - @property - def summary(self): - summary = pd.DataFrame(columns=OBJ_FIELDS) - for i, obj in enumerate(self.objectives): - for col in summary.columns: - summary.loc[i, col] = getattr(obj, col) - - # convert dtypes - for attr in ["minimize", "log"]: - summary[attr] = summary[attr].astype(bool) - - return summary - - def __repr__(self): - return self.summary.__repr__() - - # @property - # def descriptions(self) -> list: - # return [obj.description for obj in self.objectives] - - @property - def keys(self) -> list: - """ - Returns an array of the objective weights. - """ - return [obj.key for obj in self.objectives] - - @property - def weights(self) -> np.array: - """ - Returns an array of the objective weights. - """ - return np.array([obj.weight for obj in self.objectives]) - - def add(self, objective): - _validate_objectives([*self.objectives, objective]) - self.objectives.append(objective) diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py new file mode 100644 index 0000000..61c8e06 --- /dev/null +++ b/bloptools/bayesian/objectives.py @@ -0,0 +1,170 @@ +from collections.abc import Sequence +from dataclasses import dataclass +from typing import Tuple, Union + +import numpy as np +import pandas as pd + +numeric = Union[float, int] + +DEFAULT_MINIMUM_SNR = 1e2 +OBJ_FIELDS = ["description", "target", "active", "limits", "weight", "log", "n", "snr", "min_snr"] + + +class DuplicateNameError(ValueError): + ... + + +def _validate_objectives(objectives): + names = [obj.name for obj in objectives] + unique_names, counts = np.unique(names, return_counts=True) + duplicate_names = unique_names[counts > 1] + if len(duplicate_names) > 0: + raise DuplicateNameError(f"Duplicate name(s) in supplied objectives: {duplicate_names}") + + +@dataclass +class Objective: + """An objective to be used by an agent. + + Parameters + ---------- + name: str + The name of the objective. This is used as a key. + description: str + A longer description for the objective. + target: float or str + One of 'min' or 'max', or a number. The agent will respectively minimize or maximize the + objective, or target the supplied number. + log: bool + Whether to apply a log to the objective, to make it more Gaussian. + weight: float + The relative importance of this objective, used when scalarizing in multi-objective optimization. + active: bool + If True, the agent will care about this objective during optimization. + limits: tuple of floats + The range of reliable measurements for the obejctive. Outside of this, data points will be rejected. + min_snr: float + The minimum signal-to-noise ratio of the objective, used when fitting the model. + units: str + A label representing the units of the objective. + """ + + name: str + description: str = None + target: Union[float, str] = "max" + log: bool = False + weight: numeric = 1.0 + active: bool = True + limits: Tuple[numeric, numeric] = None + min_snr: numeric = DEFAULT_MINIMUM_SNR + units: str = None + + def __post_init__(self): + if self.limits is None: + if self.log: + self.limits = (0, np.inf) + else: + self.limits = (-np.inf, np.inf) + + if type(self.target) is str: + if self.target not in ["min", "max"]: + raise ValueError("'target' must be either 'min', 'max', or a number.") + + @property + def label(self): + return f"{'log ' if self.log else ''}{self.description}" + + @property + def summary(self): + series = pd.Series() + for col in OBJ_FIELDS: + series[col] = getattr(self, col) + return series + + def __repr__(self): + return self.summary.__repr__() + + @property + def noise(self): + return self.model.likelihood.noise.item() if hasattr(self, "model") else None + + @property + def snr(self): + return np.round(1 / self.model.likelihood.noise.sqrt().item(), 1) if hasattr(self, "model") else None + + @property + def n(self): + return self.model.train_targets.shape[0] if hasattr(self, "model") else 0 + + +class ObjectiveList(Sequence): + def __init__(self, objectives: list = []): + _validate_objectives(objectives) + self.objectives = objectives + + def __getitem__(self, i): + if type(i) is int: + return self.objectives[i] + elif type(i) is str: + return self.objectives[self.names.index(i)] + else: + raise ValueError(f"Invalid index {i}. An ObjectiveList must be indexed by either an integer or a string.") + + def __len__(self): + return len(self.objectives) + + @property + def summary(self): + summary = pd.DataFrame(columns=OBJ_FIELDS) + for obj in self.objectives: + for col in summary.columns: + summary.loc[obj.name, col] = getattr(obj, col) + + # convert dtypes + for attr in ["log"]: + summary[attr] = summary[attr].astype(bool) + + return summary + + def __repr__(self): + return self.summary.__repr__() + + @property + def descriptions(self) -> list: + """ + Returns an array of the objective names. + """ + return [obj.description for obj in self.objectives] + + @property + def names(self) -> list: + """ + Returns an array of the objective names. + """ + return [obj.name for obj in self.objectives] + + @property + def targets(self) -> np.array: + """ + Returns an array of the objective targets. + """ + return [obj.target for obj in self.objectives] + + @property + def weights(self) -> np.array: + """ + Returns an array of the objective weights. + """ + return np.array([obj.weight for obj in self.objectives]) + + @property + def signed_weights(self) -> np.array: + """ + Returns a signed array of the objective weights. + """ + return np.array([(1 if obj.target == "max" else -1) * obj.weight for obj in self.objectives]) + + def add(self, objective): + _validate_objectives([*self.objectives, objective]) + self.objectives.append(objective) diff --git a/bloptools/bayesian/plans.py b/bloptools/bayesian/plans.py new file mode 100644 index 0000000..26f0940 --- /dev/null +++ b/bloptools/bayesian/plans.py @@ -0,0 +1,29 @@ +import bluesky.plan_stubs as bps +import bluesky.plans as bp +import numpy as np + + +def list_scan_with_delay(*args, delay=0, **kwargs): + "Accepts all the normal 'scan' parameters, plus an optional delay." + + def one_nd_step_with_delay(detectors, step, pos_cache): + "This is a copy of bluesky.plan_stubs.one_nd_step with a sleep added." + motors = step.keys() + yield from bps.move_per_step(step, pos_cache) + yield from bps.sleep(delay) + yield from bps.trigger_and_read(list(detectors) + list(motors)) + + kwargs.setdefault("per_step", one_nd_step_with_delay) + uid = yield from bp.list_scan(*args, **kwargs) + return uid + + +def default_acquisition_plan(dofs, inputs, dets, **kwargs): + delay = kwargs.get("delay", 0) + args = [] + for dof, points in zip(dofs, np.atleast_2d(inputs).T): + args.append(dof) + args.append(list(points)) + + uid = yield from list_scan_with_delay(dets, *args, delay=delay) + return uid diff --git a/bloptools/bayesian/plotting.py b/bloptools/bayesian/plotting.py index 7e065fb..fe5fb9f 100644 --- a/bloptools/bayesian/plotting.py +++ b/bloptools/bayesian/plotting.py @@ -14,9 +14,9 @@ def _plot_objs_one_dof(agent, size=16, lw=1e0): agent.obj_fig, agent.obj_axes = plt.subplots( - agent.n_objs, + len(agent.objectives), 1, - figsize=(6, 4 * agent.n_objs), + figsize=(6, 4 * len(agent.objectives)), sharex=True, constrained_layout=True, ) @@ -27,7 +27,7 @@ def _plot_objs_one_dof(agent, size=16, lw=1e0): x_values = agent.table.loc[:, x_dof.device.name].values for obj_index, obj in enumerate(agent.objectives): - obj_fitness = agent._get_objective_targets(obj_index) + obj_values = agent.train_targets(obj.name).squeeze(-1).numpy() color = DEFAULT_COLOR_LIST[obj_index] @@ -38,7 +38,7 @@ def _plot_objs_one_dof(agent, size=16, lw=1e0): test_mean = test_posterior.mean[..., 0].detach().numpy() test_sigma = test_posterior.variance.sqrt()[..., 0].detach().numpy() - agent.obj_axes[obj_index].scatter(x_values, obj_fitness, s=size, color=color) + agent.obj_axes[obj_index].scatter(x_values, obj_values, s=size, color=color) for z in [0, 1, 2]: agent.obj_axes[obj_index].fill_between( @@ -87,12 +87,12 @@ def _plot_objs_many_dofs(agent, axes=(0, 1), shading="nearest", cmap=DEFAULT_COL test_y = test_inputs[..., 0, axes[1]].detach().squeeze().numpy() for obj_index, obj in enumerate(agent.objectives): - obj_values = agent._get_objective_targets(obj_index) + targets = agent.train_targets(obj.name).squeeze(-1).numpy() - obj_vmin, obj_vmax = np.nanpercentile(obj_values, q=[1, 99]) + obj_vmin, obj_vmax = np.nanpercentile(targets, q=[1, 99]) obj_norm = mpl.colors.Normalize(obj_vmin, obj_vmax) - data_ax = agent.obj_axes[obj_index, 0].scatter(x_values, y_values, c=obj_values, s=size, norm=obj_norm, cmap=cmap) + data_ax = agent.obj_axes[obj_index, 0].scatter(x_values, y_values, c=targets, s=size, norm=obj_norm, cmap=cmap) # mean and sigma will have shape (*input_shape,) test_posterior = obj.model.posterior(test_inputs) @@ -161,7 +161,7 @@ def _plot_objs_many_dofs(agent, axes=(0, 1), shading="nearest", cmap=DEFAULT_COL va="baseline", ) - if agent.n_objs > 1: + if len(agent.objectives) > 1: for row_index, ax in enumerate(agent.obj_axes[:, 0]): ax.annotate( agent.objectives[row_index].name, @@ -233,13 +233,14 @@ def _plot_acqf_many_dofs( # test_inputs has shape (..., 1, n_active_dofs) test_inputs = agent.test_inputs_grid() if gridded else agent.test_inputs(n=1024) + *test_dim, input_dim = test_inputs.shape test_x = test_inputs[..., 0, axes[0]].detach().squeeze().numpy() test_y = test_inputs[..., 0, axes[1]].detach().squeeze().numpy() for iacq_func, acq_func_identifier in enumerate(acq_funcs): acq_func, acq_func_meta = acquisition.get_acquisition_function(agent, acq_func_identifier) - test_acqf = acq_func(test_inputs).detach().squeeze().numpy() + test_acqf = acq_func(test_inputs.reshape(-1, 1, input_dim)).detach().reshape(test_dim).squeeze().numpy() if gridded: agent.acq_axes[iacq_func].set_title(acq_func_meta["name"]) @@ -271,7 +272,7 @@ def _plot_acqf_many_dofs( def _plot_valid_one_dof(agent, size=16, lw=1e0): - agent.valid_fig, agent.valid_ax = plt.subplots(1, 1, figsize=(6, 4 * agent.n_objs), constrained_layout=True) + agent.valid_fig, agent.valid_ax = plt.subplots(1, 1, figsize=(6, 4 * len(agent.objectives)), constrained_layout=True) x_dof = agent.dofs.subset(active=True)[0] x_values = agent.table.loc[:, x_dof.device.name].values @@ -285,7 +286,7 @@ def _plot_valid_one_dof(agent, size=16, lw=1e0): def _plot_valid_many_dofs(agent, axes=[0, 1], shading="nearest", cmap=DEFAULT_COLORMAP, size=16, gridded=None): - agent.valid_fig, agent.valid_axes = plt.subplots(1, 2, figsize=(8, 4 * agent.n_objs), constrained_layout=True) + agent.valid_fig, agent.valid_axes = plt.subplots(1, 2, figsize=(8, 4 * len(agent.objectives)), constrained_layout=True) plottable_dofs = agent.dofs.subset(active=True, read_only=False) @@ -335,9 +336,9 @@ def _plot_history(agent, x_key="index", show_all_objs=False): num_obj_plots = 1 if show_all_objs: - num_obj_plots = agent.n_objs + 1 + num_obj_plots = len(agent.objectives) + 1 - agent.n_objs + 1 if agent.n_objs > 1 else 1 + len(agent.objectives) + 1 if len(agent.objectives) > 1 else 1 hist_fig, hist_axes = plt.subplots( num_obj_plots, 1, figsize=(6, 4 * num_obj_plots), sharex=True, constrained_layout=True, dpi=200 diff --git a/bloptools/bayesian/transforms.py b/bloptools/bayesian/transforms.py new file mode 100644 index 0000000..1b28c6d --- /dev/null +++ b/bloptools/bayesian/transforms.py @@ -0,0 +1,74 @@ +from typing import Union + +import botorch +from botorch.acquisition.objective import PosteriorTransform +from botorch.posteriors.gpytorch import GPyTorchPosterior +from botorch.posteriors.posterior_list import PosteriorList +from torch import Tensor + + +class TargetingPosteriorTransform(PosteriorTransform): + r"""An affine posterior transform for scalarizing multi-output posteriors.""" + + scalarize: bool = True + + def __init__(self, weights: Tensor, targets: Tensor) -> None: + r""" + Args: + weights: A one-dimensional tensor with `m` elements representing the + linear weights on the outputs. + offset: An offset to be added to posterior mean. + """ + super().__init__() + self.targets = targets + self.register_buffer("weights", weights) + + def sampled_transform(self, y): + for i, target in enumerate(self.targets): + if target == "min": + y[..., i] = -y[..., i] + elif target != "max": + y[..., i] = -(y[..., i] - target).abs() + + return y @ self.weights.unsqueeze(-1) + + def mean_transform(self, mean, var): + for i, target in enumerate(self.targets): + if target == "min": + mean[..., i] = -mean[..., i] + elif target != "max": + mean[..., i] = -(mean[..., i] - target).abs() + + return mean @ self.weights.unsqueeze(-1) + + def variance_transform(self, mean, var): + return var @ self.weights.unsqueeze(-1) + + def evaluate(self, Y: Tensor) -> Tensor: + r"""Evaluate the transform on a set of outcomes. + + Args: + Y: A `batch_shape x q x m`-dim tensor of outcomes. + + Returns: + A `batch_shape x q`-dim tensor of transformed outcomes. + """ + return self.sampled_transform(Y) + + def forward(self, posterior: Union[GPyTorchPosterior, PosteriorList]) -> GPyTorchPosterior: + r"""Compute the posterior of the affine transformation. + + Args: + posterior: A posterior with the same number of outputs as the + elements in `self.weights`. + + Returns: + A single-output posterior. + """ + + return botorch.posteriors.transformed.TransformedPosterior( + posterior, + sample_transform=self.sampled_transform, + mean_transform=self.mean_transform, + variance_transform=self.variance_transform, + ) diff --git a/bloptools/tests/conftest.py b/bloptools/tests/conftest.py index e32b19a..c86ff9f 100644 --- a/bloptools/tests/conftest.py +++ b/bloptools/tests/conftest.py @@ -9,6 +9,7 @@ from ophyd.utils import make_dir_tree from bloptools.bayesian import DOF, Agent, Objective +from bloptools.bayesian.dofs import BrownianMotion from bloptools.utils import functions @@ -54,7 +55,7 @@ def agent(db): DOF(name="x2", limits=(-8.0, 8.0)), ] - objectives = [Objective(key="himmelblau", minimize=True)] + objectives = [Objective(name="himmelblau", target="min")] agent = Agent( dofs=dofs, @@ -71,15 +72,15 @@ def agent(db): @pytest.fixture(scope="function") def multi_agent(db): """ - A simple agent minimizing two Styblinski-Tang functions + A simple agent minimizing two Himmelblau's functions """ def digestion(db, uid): products = db[uid].table() for index, entry in products.iterrows(): - products.loc[index, "ST1"] = functions.styblinski_tang(entry.x1, entry.x2) - products.loc[index, "ST2"] = functions.styblinski_tang(entry.x1, -entry.x2) + products.loc[index, "obj1"] = functions.himmelblau(entry.x1, entry.x2) + products.loc[index, "obj2"] = functions.himmelblau(entry.x2, entry.x1) return products @@ -88,7 +89,7 @@ def digestion(db, uid): DOF(name="x2", limits=(-5.0, 5.0)), ] - objectives = [Objective(key="ST1", minimize=True), Objective(key="ST2", minimize=True)] + objectives = [Objective(name="obj1", target="min"), Objective(name="obj2", target="min")] agent = Agent( dofs=dofs, @@ -102,6 +103,36 @@ def digestion(db, uid): return agent +@pytest.fixture(scope="function") +def agent_with_passive_dofs(db): + """ + A simple agent minimizing two Himmelblau's functions + """ + + dofs = [ + DOF(name="x1", limits=(-5.0, 5.0)), + DOF(name="x2", limits=(-5.0, 5.0)), + DOF(name="x3", limits=(-5.0, 5.0), active=False), + DOF(BrownianMotion(name="brownian1"), read_only=True), + DOF(BrownianMotion(name="brownian2"), read_only=True, active=False), + ] + + objectives = [ + Objective(name="himmelblau", target="min"), + ] + + agent = Agent( + dofs=dofs, + objectives=objectives, + digestion=functions.constrained_himmelblau_digestion, + db=db, + verbose=True, + tolerate_acquisition_errors=False, + ) + + return agent + + @pytest.fixture(scope="function") def make_dirs(): root_dir = "/tmp/sirepo-bluesky-data" diff --git a/bloptools/tests/test_passive_dofs.py b/bloptools/tests/test_passive_dofs.py index 8d7ec33..797edde 100644 --- a/bloptools/tests/test_passive_dofs.py +++ b/bloptools/tests/test_passive_dofs.py @@ -1,34 +1,10 @@ import pytest -from bloptools.bayesian import DOF, Agent, BrownianMotion, Objective -from bloptools.utils import functions - @pytest.mark.test_func -def test_passive_dofs(RE, db): - dofs = [ - DOF(name="x1", limits=(-5.0, 5.0)), - DOF(name="x2", limits=(-5.0, 5.0)), - DOF(name="x3", limits=(-5.0, 5.0), active=False), - DOF(BrownianMotion(name="brownian1"), read_only=True), - DOF(BrownianMotion(name="brownian2"), read_only=True, active=False), - ] - - objectives = [ - Objective(key="himmelblau", minimize=True), - ] - - agent = Agent( - dofs=dofs, - objectives=objectives, - digestion=functions.constrained_himmelblau_digestion, - db=db, - verbose=True, - tolerate_acquisition_errors=False, - ) - - RE(agent.learn("qr", n=32)) +def test_passive_dofs(agent_with_passive_dofs, RE): + RE(agent_with_passive_dofs.learn("qr", n=32)) - agent.plot_objectives() - agent.plot_acquisition() - agent.plot_constraint() + agent_with_passive_dofs.plot_objectives() + agent_with_passive_dofs.plot_acquisition() + agent_with_passive_dofs.plot_constraint() diff --git a/bloptools/tests/test_agent.py b/bloptools/tests/test_read_write.py similarity index 89% rename from bloptools/tests/test_agent.py rename to bloptools/tests/test_read_write.py index 0c0e074..9c225d9 100644 --- a/bloptools/tests/test_agent.py +++ b/bloptools/tests/test_read_write.py @@ -9,7 +9,7 @@ def test_agent_save_load_data(agent, RE): RE(agent.learn("qr", n=16)) agent.save_data("/tmp/test_save_data.h5") agent.reset() - agent.learn(data_file="/tmp/test_save_data.h5") + agent.load_data(data_file="/tmp/test_save_data.h5") RE(agent.learn("qr", n=16)) diff --git a/bloptools/utils/functions.py b/bloptools/utils/functions.py index b46426b..810a517 100644 --- a/bloptools/utils/functions.py +++ b/bloptools/utils/functions.py @@ -114,6 +114,25 @@ def gaussian_beam_waist(x1, x2): return np.sqrt(1 + 0.25 * (x1 - x2) ** 2 + 16 * (x1 + x2 - 2) ** 2) +def hartmann4(*x): + X = np.c_[x] + + alpha = np.array([1.0, 1.2, 3.0, 3.2]) + + A = np.array([[10, 3, 17, 3.5], [0.05, 10, 17, 0.1], [3, 3.5, 1.7, 10], [17, 8, 0.05, 10]]) + + P = 1e-4 * np.array( + [ + [1312, 1696, 5569, 124], + [2329, 4135, 8307, 3736], + [2348, 1451, 3522, 2883], + [4047, 8828, 8732, 5743], + ] + ) + + return -(alpha * np.exp(-(A * np.square(X - P)).sum(axis=1))).sum() + + def hartmann6(*x): X = np.c_[x] diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/himmelblau.ipynb index f060f4a..6fb98ac 100644 --- a/docs/source/tutorials/himmelblau.ipynb +++ b/docs/source/tutorials/himmelblau.ipynb @@ -43,7 +43,7 @@ "from matplotlib import pyplot as plt\n", "from bloptools.utils import functions\n", "\n", - "x1 = x2 = np.linspace(-6, 6, 1024)\n", + "x1 = x2 = np.linspace(-6, 6, 256)\n", "X1, X2 = np.meshgrid(x1, x2)\n", "\n", "F = functions.himmelblau(X1, X2)\n", @@ -82,7 +82,7 @@ "id": "54b6f23e", "metadata": {}, "source": [ - "We also need to give the agent something to do. We want our agent to look in the feedback for a variable called \"himmelblau\", and try to minimize it." + "We also need to give the agent something to do. We want our agent to look in the feedback for a variable called 'himmelblau', and try to minimize it." ] }, { @@ -94,7 +94,7 @@ "source": [ "from bloptools.bayesian import Objective\n", "\n", - "objectives = [Objective(key=\"himmelblau\", minimize=True)]" + "objectives = [Objective(name=\"himmelblau\", description=\"Himmeblau's function\", target=\"min\")]" ] }, { @@ -246,11 +246,11 @@ }, "outputs": [], "source": [ - "X, _ = agent.ask(\"qei\", n=8, route=True)\n", + "res = agent.ask(\"qei\", n=8, route=True)\n", "agent.plot_acquisition(acq_func=\"qei\")\n", - "plt.scatter(*X.T, marker=\"d\", facecolor=\"w\", edgecolor=\"k\")\n", + "plt.scatter(*res[\"points\"].T, marker=\"d\", facecolor=\"w\", edgecolor=\"k\")\n", "plt.plot(\n", - " *X.T,\n", + " *res[\"points\"].T,\n", " color=\"r\",\n", ")" ] @@ -313,7 +313,7 @@ }, "vscode": { "interpreter": { - "hash": "eee21ccc240bdddd7cf04478199e20f7257541e2f592ca1a4d34ebdc0225d742" + "hash": "857d19a2fd370900ed798add63a0e418d98c0c9c9169a1442a8e3b86b5805755" } } }, diff --git a/docs/source/tutorials/hyperparameters.ipynb b/docs/source/tutorials/hyperparameters.ipynb index b308511..5ec7645 100644 --- a/docs/source/tutorials/hyperparameters.ipynb +++ b/docs/source/tutorials/hyperparameters.ipynb @@ -80,7 +80,7 @@ "]\n", "\n", "objectives = [\n", - " Objective(key=\"booth\", minimize=True),\n", + " Objective(name=\"booth\", target=\"min\"),\n", "]\n", "\n", "\n", @@ -145,7 +145,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.0" + "version": "3.11.5" }, "vscode": { "interpreter": { diff --git a/docs/source/tutorials/passive-dofs.ipynb b/docs/source/tutorials/passive-dofs.ipynb index 5e5a27c..71330a6 100644 --- a/docs/source/tutorials/passive-dofs.ipynb +++ b/docs/source/tutorials/passive-dofs.ipynb @@ -51,7 +51,7 @@ "]\n", "\n", "objectives = [\n", - " Objective(key=\"himmelblau\", minimize=True),\n", + " Objective(name=\"himmelblau\", target=\"min\"),\n", "]\n", "\n", "agent = Agent(\n", @@ -66,16 +66,6 @@ "RE(agent.learn(\"qr\", n=16))" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "49127601", - "metadata": {}, - "outputs": [], - "source": [ - "agent.dofs" - ] - }, { "cell_type": "code", "execution_count": null, @@ -103,7 +93,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.0" + "version": "3.11.5" }, "vscode": { "interpreter": { diff --git a/docs/wip/constrained-himmelblau copy.ipynb b/docs/wip/constrained-himmelblau copy.ipynb index 4efff8a..10b3a2c 100644 --- a/docs/wip/constrained-himmelblau copy.ipynb +++ b/docs/wip/constrained-himmelblau copy.ipynb @@ -35,7 +35,7 @@ "X1, X2 = np.meshgrid(x1, x2)\n", "from bloptools.tasks import Task\n", "\n", - "task = Task(key=\"himmelblau\", kind=\"min\")\n", + "task = Task(name=\"himmelblau\", kind=\"min\")\n", "F = functions.constrained_himmelblau(X1, X2)\n", "\n", "plt.pcolormesh(x1, x2, F, norm=mpl.colors.LogNorm(), cmap=\"gnuplot\")\n",