From 83464325fc7bb81204464d3a16ca43b0475b5c66 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 31 Oct 2023 13:20:49 -0400 Subject: [PATCH 01/17] added targeting capability to agent objectives --- bloptools/bayesian/__init__.py | 3 +- bloptools/bayesian/acquisition/analytic.py | 48 ---- bloptools/bayesian/agent.py | 16 +- bloptools/bayesian/devices.py | 212 +-------------- bloptools/bayesian/dofs.py | 249 ++++++++++++++++++ .../bayesian/{objective.py => objectives.py} | 83 ++---- bloptools/bayesian/plans.py | 29 ++ bloptools/tests/test_targeting.py | 34 +++ docs/source/tutorials/himmelblau.ipynb | 2 +- docs/source/tutorials/hyperparameters.ipynb | 2 +- docs/source/tutorials/passive-dofs.ipynb | 12 +- 11 files changed, 359 insertions(+), 331 deletions(-) create mode 100644 bloptools/bayesian/dofs.py rename bloptools/bayesian/{objective.py => objectives.py} (51%) create mode 100644 bloptools/bayesian/plans.py create mode 100644 bloptools/tests/test_targeting.py diff --git a/bloptools/bayesian/__init__.py b/bloptools/bayesian/__init__.py index 822f7be..802690e 100644 --- a/bloptools/bayesian/__init__.py +++ b/bloptools/bayesian/__init__.py @@ -1,3 +1,4 @@ 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/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..16431dc 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -22,10 +22,12 @@ 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 DOFList +from .objectives import ObjectiveList +from .plans import default_acquisition_plan + +os.environ["KMP_DUPLICATE_LIB_OK"] = "True" warnings.filterwarnings("ignore", category=botorch.exceptions.warnings.InputDataWarning) @@ -477,8 +479,14 @@ def _get_objective_targets(self, i): # transform if needed if obj.log: targets = np.where(valid, np.log(targets), np.nan) + if obj.target not in ["min", "max"]: + targets = -np.square(np.log(targets) - np.log(obj.target)) + + else: + if obj.target not in ["min", "max"]: + targets = -np.square(targets - obj.target) - if obj.minimize: + if obj.target == "min": targets *= -1 return targets diff --git a/bloptools/bayesian/devices.py b/bloptools/bayesian/devices.py index d39bc1f..3208a4d 100644 --- a/bloptools/bayesian/devices.py +++ b/bloptools/bayesian/devices.py @@ -1,218 +1,8 @@ import time as ttime import uuid -from collections.abc import Sequence -from typing import Tuple, Union import numpy as np -import pandas as pd -from ophyd import Signal, SignalRO - -DEFAULT_BOUNDS = (-5.0, +5.0) -DOF_FIELDS = ["description", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] - -numeric = Union[float, int] - - -class ReadOnlyError(Exception): - ... - - -def _validate_dofs(dofs): - """Check that a list of DOFs can be combined into a DOFList.""" - - # check that dof names are unique - unique_dof_names, counts = np.unique([dof.name for dof in dofs], 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}"') - - return list(dofs) - - -class DOF: - """A degree of freedom (DOF), to be used by an agent. - - Parameters - ---------- - name: str - The name of the DOF. This is used as a key. - description: str - A longer name 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, - 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, - ): - 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 - - 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()) - - 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 - def readback(self): - return self.device.read()[self.device.name]["value"] - - @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) - return 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 ''}" - - -class DOFList(Sequence): - """A class for handling a list of DOFs.""" - - 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 __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__() - - @property - def summary(self) -> pd.DataFrame: - table = pd.DataFrame(columns=DOF_FIELDS) - for dof in self.dofs: - for attr in table.columns: - table.loc[dof.name, attr] = getattr(dof, attr) - - # convert dtypes - for attr in ["readback", "lower_limit", "upper_limit"]: - table[attr] = table[attr].astype(float) - - for attr in ["read_only", "active"]: - table[attr] = table[attr].astype(bool) - - return table - - @property - def names(self) -> list: - return [dof.name for dof in self.dofs] - - @property - def devices(self) -> list: - return [dof.device for dof in self.dofs] - - @property - def lower_limits(self) -> np.array: - return np.array([dof.lower_limit for dof in self.dofs]) - - @property - def upper_limits(self) -> np.array: - return np.array([dof.upper_limit for dof in self.dofs]) - - @property - def limits(self) -> np.array: - """ - Returns a (n_dof, 2) array of bounds. - """ - return np.c_[self.lower_limits, self.upper_limits] - - @property - def readback(self) -> np.array: - return np.array([dof.readback for dof in self.dofs]) - - def add(self, dof): - _validate_dofs([*self.dofs, dof]) - self.dofs.append(dof) - - def _dof_read_only_mask(self, read_only=None): - return [dof.read_only == read_only if read_only is not None else True for dof in self.dofs] - - def _dof_active_mask(self, active=None): - return [dof.active == active if active is not None else True for dof in self.dofs] - - def _dof_tags_mask(self, tags=[]): - return [np.isin(dof["tags"], tags).any() if tags else True for dof in self.dofs] - - def _dof_mask(self, active=None, read_only=None, tags=[]): - return [ - (k and m and t) - for k, m, t in zip(self._dof_read_only_mask(read_only), self._dof_active_mask(active), self._dof_tags_mask(tags)) - ] - - 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 +from ophyd import Signal, SignalRO # noqa F401 class BrownianMotion(SignalRO): diff --git a/bloptools/bayesian/dofs.py b/bloptools/bayesian/dofs.py new file mode 100644 index 0000000..33eda63 --- /dev/null +++ b/bloptools/bayesian/dofs.py @@ -0,0 +1,249 @@ +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 +import pandas as pd +from ophyd import Signal, SignalRO + +DEFAULT_BOUNDS = (-5.0, +5.0) +DOF_FIELDS = ["name", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] + +numeric = Union[float, int] + + +class ReadOnlyError(Exception): + ... + + +def _validate_dofs(dofs): + dof_names = [dof.name for dof in dofs] + + # check that dof names are unique + 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}"') + + return list(dofs) + + +@dataclass +class DOF: + device: Signal = None + limits: Tuple[float, float] = (-10.0, 10.0) + name: str = None + units: str = "" + read_only: bool = False + active: bool = True + tags: list = field(default_factory=list) + latent_group = None + + def __post_init__(self): + self.uuid = str(uuid.uuid4()) + + if self.name is None: + self.name = self.device.name if hasattr(self.device, "name") else self.uuid + + 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): + return float(self.limits[0]) + + @property + def upper_limit(self): + return float(self.limits[1]) + + @property + def readback(self): + return self.device.read()[self.device.name]["value"] + + @property + def summary(self) -> pd.Series: + series = pd.Series(index=DOF_FIELDS) + for attr in series.index: + series[attr] = getattr(self, attr) + return series + + @property + def label(self) -> str: + 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): + def __init__(self, dofs: list = []): + _validate_dofs(dofs) + self.dofs = dofs + + def __getitem__(self, i): + return self.dofs[i] + + def __len__(self): + return len(self.dofs) + + def __repr__(self): + 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) + for i, dof in enumerate(self.dofs): + for attr in table.columns: + table.loc[i, attr] = getattr(dof, attr) + + # convert dtypes + for attr in ["readback", "lower_limit", "upper_limit"]: + table[attr] = table[attr].astype(float) + + for attr in ["read_only", "active"]: + table[attr] = table[attr].astype(bool) + + return table + + @property + def names(self) -> list: + return [dof.name for dof in self.dofs] + + @property + 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]) + + @property + def upper_limits(self) -> np.array: + return np.array([dof.upper_limit for dof in self.dofs]) + + @property + def limits(self) -> np.array: + """ + Returns a (n_dof, 2) array of bounds. + """ + return np.c_[self.lower_limits, self.upper_limits] + + @property + def readback(self) -> np.array: + return np.array([dof.readback for dof in self.dofs]) + + def add(self, dof): + _validate_dofs([*self.dofs, dof]) + self.dofs.append(dof) + + def _dof_read_only_mask(self, read_only=None): + return [dof.read_only == read_only if read_only is not None else True for dof in self.dofs] + + def _dof_active_mask(self, active=None): + return [dof.active == active if active is not None else True for dof in self.dofs] + + def _dof_tags_mask(self, tags=[]): + return [np.isin(dof["tags"], tags).any() if tags else True for dof in self.dofs] + + def _dof_mask(self, active=None, read_only=None, tags=[]): + return [ + (k and m and t) + for k, m, t in zip(self._dof_read_only_mask(read_only), self._dof_active_mask(active), self._dof_tags_mask(tags)) + ] + + 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 _subset_devices(self, read_only=None, active=None, tags=[]): + # return [dof["device"] for dof in self._subset_dofs(read_only, active, tags)] + + # def _read_subset_devices(self, read_only=None, active=None, tags=[]): + # return [device.read()[device.name]["value"] for device in self._subset_devices(read_only, active, tags)] + + # def _subset_dof_names(self, read_only=None, active=None, tags=[]): + # return [device.name for device in self._subset_devices(read_only, active, tags)] + + # def _subset_dof_limits(self, read_only=None, active=None, tags=[]): + # dofs_subset = self._subset_dofs(read_only, active, tags) + # if len(dofs_subset) > 0: + # return torch.tensor([dof["limits"] for dof in dofs_subset], dtype=torch.float64).T + # return torch.empty((2, 0)) + + def activate(self, read_only=None, active=None, tags=[]): + for dof in self._subset_dofs(read_only, active, tags): + dof.active = True + + def deactivate(self, read_only=None, active=None, tags=[]): + for dof in self._subset_dofs(read_only, active, tags): + dof.active = False + + +class BrownianMotion(SignalRO): + """ + Read-only degree of freedom simulating brownian motion + """ + + def __init__(self, name=None, theta=0.95, *args, **kwargs): + name = name if name is not None else str(uuid.uuid4()) + + super().__init__(name=name, *args, **kwargs) + + self.theta = theta + self.old_t = ttime.monotonic() + self.old_y = 0.0 + + def get(self): + new_t = ttime.monotonic() + alpha = self.theta ** (new_t - self.old_t) + new_y = alpha * self.old_y + np.sqrt(1 - alpha**2) * np.random.standard_normal() + + self.old_t = new_t + self.old_y = new_y + return new_y + + +class TimeReadback(SignalRO): + """ + Returns the current timestamp. + """ + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def get(self): + return ttime.time() + + +class ConstantReadback(SignalRO): + """ + Returns a constant every time you read it (more useful than you'd think). + """ + + def __init__(self, constant=1, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.constant = constant + + def get(self): + return self.constant diff --git a/bloptools/bayesian/objective.py b/bloptools/bayesian/objectives.py similarity index 51% rename from bloptools/bayesian/objective.py rename to bloptools/bayesian/objectives.py index 47b963f..9944a4b 100644 --- a/bloptools/bayesian/objective.py +++ b/bloptools/bayesian/objectives.py @@ -1,8 +1,10 @@ from collections.abc import Sequence +from dataclasses import dataclass from typing import Tuple, Union import numpy as np import pandas as pd +from ophyd import Signal numeric = Union[float, int] @@ -22,63 +24,36 @@ def _validate_objectives(objectives): raise DuplicateKeyError(f'Duplicate key(s) in supplied objectives: "{duplicate_keys}"') +@dataclass 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) + key: str + name: str = None + target: float | str = "max" + log: bool = False + weight: numeric = 1.0 + limits: Tuple[numeric, numeric] = None + min_snr: numeric = DEFAULT_MINIMUM_SNR + units: str = None + + def __post_init__(self): + if self.name is None: + self.name = self.key + + 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.") + + self.device = Signal(name=self.name) @property def label(self): - return f"{'neg ' if self.minimize else ''}{'log ' if self.log else ''}{self.description}" + return f"{'neg ' if self.target == 'min' else ''}{'log ' if self.log else ''}{self.name}" @property def summary(self): @@ -114,7 +89,7 @@ def summary(self): summary.loc[i, col] = getattr(obj, col) # convert dtypes - for attr in ["minimize", "log"]: + for attr in ["log"]: summary[attr] = summary[attr].astype(bool) return summary 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/tests/test_targeting.py b/bloptools/tests/test_targeting.py new file mode 100644 index 0000000..3db1c06 --- /dev/null +++ b/bloptools/tests/test_targeting.py @@ -0,0 +1,34 @@ +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)) + + agent.plot_objectives() + agent.plot_acquisition() + agent.plot_validity() diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/himmelblau.ipynb index f060f4a..1ac5c6b 100644 --- a/docs/source/tutorials/himmelblau.ipynb +++ b/docs/source/tutorials/himmelblau.ipynb @@ -94,7 +94,7 @@ "source": [ "from bloptools.bayesian import Objective\n", "\n", - "objectives = [Objective(key=\"himmelblau\", minimize=True)]" + "objectives = [Objective(key=\"himmelblau\", target=\"min\")]" ] }, { diff --git a/docs/source/tutorials/hyperparameters.ipynb b/docs/source/tutorials/hyperparameters.ipynb index b308511..02f2a9f 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(key=\"booth\", target=\"min\"),\n", "]\n", "\n", "\n", diff --git a/docs/source/tutorials/passive-dofs.ipynb b/docs/source/tutorials/passive-dofs.ipynb index 5e5a27c..3af9c40 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(key=\"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, From 380ee116fa0456529197bdfb1be19b6501787b65 Mon Sep 17 00:00:00 2001 From: XPD Operator Date: Tue, 31 Oct 2023 15:18:16 -0400 Subject: [PATCH 02/17] work at xpd --- bloptools/bayesian/agent.py | 43 ++++++++------------------------ bloptools/bayesian/objectives.py | 2 +- 2 files changed, 12 insertions(+), 33 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 16431dc..f692da6 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -139,7 +139,7 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= 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) + skew_dims = self.latent_dim_tuples 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 @@ -234,7 +234,7 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): print(f'finding points with acquisition function "{acq_func_name}" ...') if acq_func_type in ["analytic", "monte_carlo"]: - if not self.initialized: + if not self.has_models: raise RuntimeError( f'Can\'t construct non-trivial acquisition function "{acq_func_identifier}"' f" (the agent is not initialized!)" @@ -361,40 +361,18 @@ def learn( 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)) - - Parameters - ---------- - acq_func : str - A valid identifier for an implemented acquisition function. - n : int - How many points to sample on each iteration. - iterations: int - How many iterations of the learning loop to perform. - train: bool - Whether to train the models upon telling the agent. - append: bool - If `True`, add the new data to the old data. If `False`, replace the old data with the new data. - data_file: str - If supplied, read a saved data file instead of running the acquisition plan. - hypers_file: str - If supplied, read a saved hyperparameter file instead of fitting models. NOTE: The agent will assume these - hyperparameters a priori for the rest of the run, and not try to fit a model. + """ + This iterates the learning algorithm, looping over ask -> acquire -> tell. + It should be passed to a Bluesky RunEngine. """ if data_file is not None: new_table = pd.read_hdf(data_file, key="table") - 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" + if self.sample_center_on_init and not self.has_models: + new_table = yield from self.acquire(self.dofs.subset(active=True, read_only=False).limits.mean(axis=1)) + new_table.loc[:, "acq_func"] = "sample_center_on_init" + self.tell(new_table=new_table, train=False) for i in range(iterations): print(f"running iteration {i + 1} / {iterations}") @@ -423,7 +401,8 @@ 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: + del obj.model def benchmark( self, output_dir="./", runs=16, n_init=64, learning_kwargs_list=[{"acq_func": "qei", "n": 4, "iterations": 16}] diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 9944a4b..6c63af8 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -28,7 +28,7 @@ def _validate_objectives(objectives): class Objective: key: str name: str = None - target: float | str = "max" + target: Union[float, str] = "max" log: bool = False weight: numeric = 1.0 limits: Tuple[numeric, numeric] = None From ef03e1487327b49672f1a8d48f61d77310079156 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 31 Oct 2023 15:29:40 -0400 Subject: [PATCH 03/17] update syntax in tests --- bloptools/tests/conftest.py | 4 ++-- bloptools/tests/test_passive_dofs.py | 2 +- bloptools/tests/test_targeting.py | 34 ---------------------------- 3 files changed, 3 insertions(+), 37 deletions(-) delete mode 100644 bloptools/tests/test_targeting.py diff --git a/bloptools/tests/conftest.py b/bloptools/tests/conftest.py index e32b19a..fa88e09 100644 --- a/bloptools/tests/conftest.py +++ b/bloptools/tests/conftest.py @@ -54,7 +54,7 @@ def agent(db): DOF(name="x2", limits=(-8.0, 8.0)), ] - objectives = [Objective(key="himmelblau", minimize=True)] + objectives = [Objective(key="himmelblau", target="min")] agent = Agent( dofs=dofs, @@ -88,7 +88,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(key="ST1", target="min"), Objective(key="ST2", target="min")] agent = Agent( dofs=dofs, diff --git a/bloptools/tests/test_passive_dofs.py b/bloptools/tests/test_passive_dofs.py index 8d7ec33..5b13b87 100644 --- a/bloptools/tests/test_passive_dofs.py +++ b/bloptools/tests/test_passive_dofs.py @@ -15,7 +15,7 @@ def test_passive_dofs(RE, db): ] objectives = [ - Objective(key="himmelblau", minimize=True), + Objective(key="himmelblau", target="min"), ] agent = Agent( diff --git a/bloptools/tests/test_targeting.py b/bloptools/tests/test_targeting.py deleted file mode 100644 index 3db1c06..0000000 --- a/bloptools/tests/test_targeting.py +++ /dev/null @@ -1,34 +0,0 @@ -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)) - - agent.plot_objectives() - agent.plot_acquisition() - agent.plot_validity() From cc0e8854f350a6274df7262ab3a47ec75f0cc49e Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Wed, 8 Nov 2023 10:24:20 -0800 Subject: [PATCH 04/17] fix bug in learn method --- bloptools/bayesian/agent.py | 84 +++++++++++++++++++------------- bloptools/bayesian/objectives.py | 10 +++- bloptools/tests/test_agent.py | 2 +- 3 files changed, 61 insertions(+), 35 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index f692da6..508374c 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -23,12 +23,10 @@ from .. import utils from . import acquisition, models, plotting from .digestion import default_digestion_function -from .dofs import DOFList -from .objectives import ObjectiveList +from .dofs import DOF, DOFList +from .objectives import Objective, ObjectiveList from .plans import default_acquisition_plan -os.environ["KMP_DUPLICATE_LIB_OK"] = "True" - warnings.filterwarnings("ignore", category=botorch.exceptions.warnings.InputDataWarning) mpl.rc("image", cmap="coolwarm") @@ -139,7 +137,9 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= self.table = pd.concat([self.table, new_table]) if append else new_table self.table.index = np.arange(len(self.table)) - skew_dims = self.latent_dim_tuples + # TODO: should be a check per model + if len(self.table) > 2: + self._update_models(train=train_models, a_priori_hypers=hypers) 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 @@ -234,7 +234,7 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): print(f'finding points with acquisition function "{acq_func_name}" ...') if acq_func_type in ["analytic", "monte_carlo"]: - if not self.has_models: + 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" (the agent is not initialized!)" @@ -350,47 +350,65 @@ def acquire(self, acquisition_inputs): return products + def load_data(self, data_file, append=True, train_models=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.keys} + metadata = new_table.to_dict(orient="list") + self.tell(x=x, y=y, metadata=metadata, append=append, train_models=train_models) + def learn( self, acq_func=None, n=1, iterations=1, upsample=1, - train=True, - data_file=None, + train_models=True, hypers_file=None, append=True, ): - """ - This iterates the learning algorithm, looping over ask -> acquire -> tell. - It should be passed to a Bluesky RunEngine. - """ - - if data_file is not None: - new_table = pd.read_hdf(data_file, key="table") - - if self.sample_center_on_init and not self.has_models: - new_table = yield from self.acquire(self.dofs.subset(active=True, read_only=False).limits.mean(axis=1)) - new_table.loc[:, "acq_func"] = "sample_center_on_init" - self.tell(new_table=new_table, train=False) + """This returns a Bluesky plan which iterates the learning algorithm, looping over ask -> acquire -> tell. - 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"] + For example: - else: - raise ValueError("You must supply either an acquisition function or a filepath!") + RE(agent.learn("qr", n=16)) + RE(agent.learn("qei", n=4, iterations=4)) - 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") + Parameters + ---------- + acq_func : str + A valid identifier for an implemented acquisition function. + n : int + How many points to sample on each iteration. + iterations: int + How many iterations of the learning loop to perform. + train: bool + Whether to train the models upon telling the agent. + append: bool + If `True`, add the new data to the old data. If `False`, replace the old data with the new data. + data_file: str + If supplied, read a saved data file instead of running the acquisition plan. + hypers_file: str + If supplied, read a saved hyperparameter file instead of fitting models. NOTE: The agent will assume these + hyperparameters a priori for the rest of the run, and not try to fit a model. + """ - self.tell(x=x, y=y, metadata=metadata, append=append, train_models=True) + 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" - self.initialized = True + for i in range(iterations): + print(f"running iteration {i + 1} / {iterations}") + for single_acq_func in np.atleast_1d(acq_func): + acq_points, acq_func_meta = self.ask(n=n, acq_func_identifier=single_acq_func) + new_table = yield from self.acquire(acq_points) + new_table.loc[:, "acq_func"] = acq_func_meta["name"] + + 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=train_models) def get_acquisition_function(self, identifier, return_metadata=False): """Returns a BoTorch acquisition function for a given identifier. Acquisition functions can be diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 6c63af8..55e4f8d 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -9,7 +9,7 @@ numeric = Union[float, int] DEFAULT_MINIMUM_SNR = 1e1 -OBJ_FIELDS = ["description", "key", "limits", "weight", "minimize", "log"] +OBJ_FIELDS = ["description", "key", "limits", "weight", "minimize", "log", "n", "snr", "min_snr"] class DuplicateKeyError(ValueError): @@ -69,6 +69,14 @@ def __repr__(self): 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 = []): diff --git a/bloptools/tests/test_agent.py b/bloptools/tests/test_agent.py index 0c0e074..9c225d9 100644 --- a/bloptools/tests/test_agent.py +++ b/bloptools/tests/test_agent.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)) From e3e18e53d3f23d13ca52b8e5ef745e12803a115c Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Wed, 8 Nov 2023 15:46:39 -0800 Subject: [PATCH 05/17] added targeting posterior transform --- bloptools/bayesian/acquisition/__init__.py | 13 ++- bloptools/bayesian/agent.py | 105 ++++++++++++++------- bloptools/bayesian/dofs.py | 7 +- bloptools/bayesian/models.py | 16 ++++ bloptools/bayesian/objectives.py | 44 +++++---- bloptools/bayesian/plotting.py | 22 ++--- bloptools/bayesian/transforms.py | 57 +++++++++++ 7 files changed, 194 insertions(+), 70 deletions(-) create mode 100644 bloptools/bayesian/transforms.py diff --git a/bloptools/bayesian/acquisition/__init__.py b/bloptools/bayesian/acquisition/__init__.py index 854e63f..ac3c1b5 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 @@ -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/agent.py b/bloptools/bayesian/agent.py index 508374c..f3c9c01 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -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 @@ -26,6 +27,7 @@ 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) @@ -144,20 +146,18 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= 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 - # if self.initialized: - # cached_hypers = self.hypers - inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) 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) + values = self.get_objective_targets(i) + + 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_targets = torch.tensor(targets[train_index], dtype=torch.double).unsqueeze(-1) # .unsqueeze(0) + train_values = torch.tensor(values[train_index], dtype=torch.double).unsqueeze(-1) # .unsqueeze(0) # for constructing the log normal noise prior # target_snr = 2e2 @@ -176,7 +176,7 @@ def _update_models(self, train=True, skew_dims=None, a_priori_hypers=None): obj.model = models.LatentGP( train_inputs=train_inputs, - train_targets=train_targets, + train_targets=train_values, likelihood=likelihood, skew_dims=skew_dims, input_transform=self._subset_input_transform(active=True), @@ -343,7 +343,7 @@ 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!") @@ -353,7 +353,7 @@ def acquire(self, acquisition_inputs): def load_data(self, data_file, append=True, train_models=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.keys} + 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_models=train_models) @@ -406,7 +406,7 @@ def learn( new_table.loc[:, "acq_func"] = acq_func_meta["name"] 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} + 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_models=train_models) @@ -459,15 +459,20 @@ 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.""" + def get_objective_targets(self, i): + """Returns the values associated with each objective.""" + obj = self.objectives[i] - targets = self.table.loc[:, obj.key].values.copy() + 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]) @@ -475,33 +480,69 @@ def _get_objective_targets(self, i): # transform if needed if obj.log: - targets = np.where(valid, np.log(targets), np.nan) - if obj.target not in ["min", "max"]: - targets = -np.square(np.log(targets) - np.log(obj.target)) + targets = np.where(targets > 0, np.log(targets), np.nan) - else: - if obj.target not in ["min", "max"]: - targets = -np.square(targets - obj.target) + return targets - if obj.target == "min": - targets *= -1 + # def _get_objective_targets(self, i): + # """Returns the targets (what we fit to) for an objective, given the objective index.""" + # obj = self.objectives[i] - return targets + # 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(valid, np.log(targets), np.nan) + # if obj.target not in ["min", "max"]: + # targets = -np.square(np.log(targets) - np.log(obj.target)) + + # else: + # if obj.target not in ["min", "max"]: + # targets = -np.square(targets - obj.target) + + # if obj.target == "min": + # targets *= -1 + + # return targets + + @property + def scalarizing_transform(self): + return ScalarizedPosteriorTransform(weights=self.objective_weights_torch, offset=0) + + @property + def targeting_transform(self): + return TargetingPosteriorTransform(weights=self.objective_weights_torch, targets=self.pseudo_targets) @property - def n_objs(self): - """Returns a (num_objectives x n_observations) array of objectives""" - return len(self.objectives) + def pseudo_targets(self): + """Targets for the posterior transform""" + return torch.tensor( + [ + self.objectives_targets[..., i].max() + if t == "max" + else self.objectives_targets[..., i].min() + if t == "min" + else t + for i, t in enumerate(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) + return torch.cat( + [torch.tensor(self.get_objective_targets(i))[..., None] for i in range(len(self.objectives))], dim=1 + ) @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.objectives_targets).sum(axis=-1) + # return (self.objectives_targets * self.objectives.signed_weights).sum(axis=-1) @property def max_scalarized_objective(self): @@ -615,7 +656,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 @@ -625,9 +666,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 @@ -635,7 +676,7 @@ def save_hypers(self, filepath): """Save the agent's fitted hyperparameters to a given filepath.""" hypers = self.hypers with h5py.File(filepath, "w") as f: - for model_key in hypers.keys(): + for model_key in hypers.names(): f.create_group(model_key) for param_key, param_value in hypers[model_key].items(): f[model_key].create_dataset(param_key, data=param_value) @@ -645,7 +686,7 @@ def load_hypers(filepath): """Load hyperparameters from a file.""" hypers = {} with h5py.File(filepath, "r") as f: - for model_key in f.keys(): + for model_key in f.names(): hypers[model_key] = OrderedDict() for param_key, param_value in f[model_key].items(): hypers[model_key][param_key] = torch.tensor(np.atleast_1d(param_value[()])) diff --git a/bloptools/bayesian/dofs.py b/bloptools/bayesian/dofs.py index 33eda63..a68651f 100644 --- a/bloptools/bayesian/dofs.py +++ b/bloptools/bayesian/dofs.py @@ -9,7 +9,7 @@ from ophyd import Signal, SignalRO DEFAULT_BOUNDS = (-5.0, +5.0) -DOF_FIELDS = ["name", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] +DOF_FIELDS = ["name", "description", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] numeric = Union[float, int] @@ -33,13 +33,14 @@ def _validate_dofs(dofs): @dataclass class DOF: device: Signal = None - limits: Tuple[float, float] = (-10.0, 10.0) + 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) - latent_group = None + latent_group: str = None def __post_init__(self): self.uuid = str(uuid.uuid4()) diff --git a/bloptools/bayesian/models.py b/bloptools/bayesian/models.py index eff9a64..fab5f38 100644 --- a/bloptools/bayesian/models.py +++ b/bloptools/bayesian/models.py @@ -21,6 +21,22 @@ def __init__(self, train_inputs, train_targets, skew_dims=True, *args, **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, + diag_prior=True, + scale=True, + **kwargs + ) + + class LatentDirichletClassifier(LatentGP): def __init__(self, train_inputs, train_targets, skew_dims=True, *args, **kwargs): super().__init__(train_inputs, train_targets, skew_dims, *args, **kwargs) diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 55e4f8d..3e13978 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -4,30 +4,29 @@ import numpy as np import pandas as pd -from ophyd import Signal numeric = Union[float, int] DEFAULT_MINIMUM_SNR = 1e1 -OBJ_FIELDS = ["description", "key", "limits", "weight", "minimize", "log", "n", "snr", "min_snr"] +OBJ_FIELDS = ["description", "limits", "weight", "minimize", "log", "n", "snr", "min_snr"] -class DuplicateKeyError(ValueError): +class DuplicateNameError(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}"') + 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: - key: str - name: str = None + name: str + description: str = None target: Union[float, str] = "max" log: bool = False weight: numeric = 1.0 @@ -36,9 +35,6 @@ class Objective: units: str = None def __post_init__(self): - if self.name is None: - self.name = self.key - if self.limits is None: if self.log: self.limits = (0, np.inf) @@ -49,11 +45,11 @@ def __post_init__(self): if self.target not in ["min", "max"]: raise ValueError("'target' must be either 'min', 'max', or a number.") - self.device = Signal(name=self.name) + # self.device = Signal(name=self.name) @property def label(self): - return f"{'neg ' if self.target == 'min' else ''}{'log ' if self.log else ''}{self.name}" + return f"{'log ' if self.log else ''}{self.description}" @property def summary(self): @@ -110,11 +106,18 @@ def __repr__(self): # return [obj.description for obj in self.objectives] @property - def keys(self) -> list: + def names(self) -> list: + """ + Returns an array of the objective weights. + """ + return [obj.name for obj in self.objectives] + + @property + def targets(self) -> np.array: """ Returns an array of the objective weights. """ - return [obj.key for obj in self.objectives] + return [obj.target for obj in self.objectives] @property def weights(self) -> np.array: @@ -123,6 +126,13 @@ def weights(self) -> np.array: """ 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/plotting.py b/bloptools/bayesian/plotting.py index 7e065fb..e774c95 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_fitness = agent.get_objective_targets(obj_index) color = DEFAULT_COLOR_LIST[obj_index] @@ -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.get_objective_targets(obj_index) - 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, @@ -271,7 +271,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 +285,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 +335,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..2f02ff8 --- /dev/null +++ b/bloptools/bayesian/transforms.py @@ -0,0 +1,57 @@ +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 # pragma: no cover +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.register_buffer("targets", targets) + self.register_buffer("weights", weights) + + self.sampled_transform = lambda y: -(y - self.targets).abs() @ self.weights.unsqueeze(-1) + self.mean_transform = lambda mean, var: -(mean - self.targets).abs() @ self.weights.unsqueeze(-1) + self.variance_transform = lambda mean, var: -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, + ) From 440bc163d1d6860e0edd358555e86788cb931373 Mon Sep 17 00:00:00 2001 From: Tanny Date: Thu, 9 Nov 2023 11:08:17 -0800 Subject: [PATCH 06/17] als tweaks --- bloptools/bayesian/agent.py | 13 +++++++------ bloptools/bayesian/objectives.py | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index f3c9c01..24f0aba 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -522,12 +522,13 @@ def pseudo_targets(self): """Targets for the posterior transform""" return torch.tensor( [ - self.objectives_targets[..., i].max() - if t == "max" - else self.objectives_targets[..., i].min() - if t == "min" - else t - for i, t in enumerate(self.objectives.targets) + 1.e32 + if obj.target == "max" + else -1.e32 + if obj.target == "min" + else np.log(obj.target) if obj.log + else obj.target + for i, obj in enumerate(self.objectives) ] ) diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 3e13978..4eeb8da 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -8,7 +8,7 @@ numeric = Union[float, int] DEFAULT_MINIMUM_SNR = 1e1 -OBJ_FIELDS = ["description", "limits", "weight", "minimize", "log", "n", "snr", "min_snr"] +OBJ_FIELDS = ["description", "target", "limits", "weight", "log", "n", "snr", "min_snr"] class DuplicateNameError(ValueError): From 6576bc81e2b0c6fd97bcbba37bb747e9dccd8f12 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Fri, 10 Nov 2023 09:55:52 -0800 Subject: [PATCH 07/17] better targeting method --- bloptools/bayesian/agent.py | 60 ++++++++------------- bloptools/bayesian/objectives.py | 4 +- bloptools/bayesian/transforms.py | 25 +++++++-- bloptools/tests/conftest.py | 8 +-- bloptools/tests/test_passive_dofs.py | 2 +- docs/source/tutorials/himmelblau.ipynb | 6 +-- docs/source/tutorials/hyperparameters.ipynb | 4 +- docs/source/tutorials/passive-dofs.ipynb | 4 +- docs/wip/constrained-himmelblau copy.ipynb | 2 +- 9 files changed, 57 insertions(+), 58 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 24f0aba..18f4215 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -49,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. @@ -107,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 @@ -135,12 +137,17 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= A dict of hyperparameters for the model to assume a priori. """ + # n_before_tell = len(self.table) + new_table = pd.DataFrame({**x, **y, **metadata} if metadata is not None else {**x, **y}) self.table = pd.concat([self.table, new_table]) if append else new_table self.table.index = np.arange(len(self.table)) + # n_after_tell = len(self.table) + # TODO: should be a check per model if len(self.table) > 2: + # if n_before_tell % self.train_every != n_after_tell % self.train_every: self._update_models(train=train_models, a_priori_hypers=hypers) def _update_models(self, train=True, skew_dims=None, a_priori_hypers=None): @@ -257,13 +264,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 @@ -484,51 +492,25 @@ def get_objective_targets(self, i): return targets - # 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.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(valid, np.log(targets), np.nan) - # if obj.target not in ["min", "max"]: - # targets = -np.square(np.log(targets) - np.log(obj.target)) - - # else: - # if obj.target not in ["min", "max"]: - # targets = -np.square(targets - obj.target) - - # if obj.target == "min": - # targets *= -1 - - # return targets - @property def scalarizing_transform(self): return ScalarizedPosteriorTransform(weights=self.objective_weights_torch, offset=0) @property def targeting_transform(self): - return TargetingPosteriorTransform(weights=self.objective_weights_torch, targets=self.pseudo_targets) + return TargetingPosteriorTransform(weights=self.objective_weights_torch, targets=self.objectives.targets) @property def pseudo_targets(self): """Targets for the posterior transform""" return torch.tensor( [ - 1.e32 - if obj.target == "max" - else -1.e32 - if obj.target == "min" - else np.log(obj.target) if obj.log - else obj.target - for i, obj in enumerate(self.objectives) + self.objectives_targets[..., i].max() + if t == "max" + else self.objectives_targets[..., i].min() + if t == "min" + else t + for i, t in enumerate(self.objectives.targets) ] ) @@ -677,7 +659,7 @@ def save_hypers(self, filepath): """Save the agent's fitted hyperparameters to a given filepath.""" hypers = self.hypers with h5py.File(filepath, "w") as f: - for model_key in hypers.names(): + for model_key in hypers.keys(): f.create_group(model_key) for param_key, param_value in hypers[model_key].items(): f[model_key].create_dataset(param_key, data=param_value) @@ -687,7 +669,7 @@ def load_hypers(filepath): """Load hyperparameters from a file.""" hypers = {} with h5py.File(filepath, "r") as f: - for model_key in f.names(): + for model_key in f.keys(): hypers[model_key] = OrderedDict() for param_key, param_value in f[model_key].items(): hypers[model_key][param_key] = torch.tensor(np.atleast_1d(param_value[()])) diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 4eeb8da..7f15752 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -88,9 +88,9 @@ def __len__(self): @property def summary(self): summary = pd.DataFrame(columns=OBJ_FIELDS) - for i, obj in enumerate(self.objectives): + for obj in self.objectives: for col in summary.columns: - summary.loc[i, col] = getattr(obj, col) + summary.loc[obj.name, col] = getattr(obj, col) # convert dtypes for attr in ["log"]: diff --git a/bloptools/bayesian/transforms.py b/bloptools/bayesian/transforms.py index 2f02ff8..9f64d9c 100644 --- a/bloptools/bayesian/transforms.py +++ b/bloptools/bayesian/transforms.py @@ -20,12 +20,29 @@ def __init__(self, weights: Tensor, targets: Tensor) -> None: offset: An offset to be added to posterior mean. """ super().__init__() - self.register_buffer("targets", targets) + self.targets = targets self.register_buffer("weights", weights) - self.sampled_transform = lambda y: -(y - self.targets).abs() @ self.weights.unsqueeze(-1) - self.mean_transform = lambda mean, var: -(mean - self.targets).abs() @ self.weights.unsqueeze(-1) - self.variance_transform = lambda mean, var: -var @ self.weights.unsqueeze(-1) + 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. diff --git a/bloptools/tests/conftest.py b/bloptools/tests/conftest.py index fa88e09..0fad97d 100644 --- a/bloptools/tests/conftest.py +++ b/bloptools/tests/conftest.py @@ -54,7 +54,7 @@ def agent(db): DOF(name="x2", limits=(-8.0, 8.0)), ] - objectives = [Objective(key="himmelblau", target="min")] + objectives = [Objective(name="himmelblau", target="min")] agent = Agent( dofs=dofs, @@ -78,8 +78,8 @@ 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 +88,7 @@ def digestion(db, uid): DOF(name="x2", limits=(-5.0, 5.0)), ] - objectives = [Objective(key="ST1", target="min"), Objective(key="ST2", target="min")] + objectives = [Objective(name="obj1", target="min"), Objective(name="obj2", target="min")] agent = Agent( dofs=dofs, diff --git a/bloptools/tests/test_passive_dofs.py b/bloptools/tests/test_passive_dofs.py index 5b13b87..92cfda4 100644 --- a/bloptools/tests/test_passive_dofs.py +++ b/bloptools/tests/test_passive_dofs.py @@ -15,7 +15,7 @@ def test_passive_dofs(RE, db): ] objectives = [ - Objective(key="himmelblau", target="min"), + Objective(name="himmelblau", target="min"), ] agent = Agent( diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/himmelblau.ipynb index 1ac5c6b..29ceb16 100644 --- a/docs/source/tutorials/himmelblau.ipynb +++ b/docs/source/tutorials/himmelblau.ipynb @@ -94,7 +94,7 @@ "source": [ "from bloptools.bayesian import Objective\n", "\n", - "objectives = [Objective(key=\"himmelblau\", target=\"min\")]" + "objectives = [Objective(name=\"himmelblau\", target=\"min\")]" ] }, { @@ -295,7 +295,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3.10.0 ('nsls2')", "language": "python", "name": "python3" }, @@ -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 02f2a9f..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\", target=\"min\"),\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 3af9c40..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\", target=\"min\"),\n", + " Objective(name=\"himmelblau\", target=\"min\"),\n", "]\n", "\n", "agent = Agent(\n", @@ -93,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", From 34323e719e93e31f429d156d77d247d8e34ff80d Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Fri, 10 Nov 2023 14:50:07 -0800 Subject: [PATCH 08/17] make objectives toggleable --- bloptools/bayesian/acquisition/__init__.py | 2 +- bloptools/bayesian/agent.py | 20 +++++++++++--------- bloptools/bayesian/dofs.py | 6 +++--- bloptools/bayesian/objectives.py | 5 +++-- bloptools/bayesian/plotting.py | 3 ++- bloptools/tests/conftest.py | 2 +- 6 files changed, 21 insertions(+), 17 deletions(-) diff --git a/bloptools/bayesian/acquisition/__init__.py b/bloptools/bayesian/acquisition/__init__.py index ac3c1b5..a66afbb 100644 --- a/bloptools/bayesian/acquisition/__init__.py +++ b/bloptools/bayesian/acquisition/__init__.py @@ -37,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 diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 18f4215..8c6dbb9 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -117,6 +117,8 @@ def __init__( self.initialized = False self.a_priori_hypers = None + self.n_last_trained = 0 + def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models=True, hypers=None): """ Inform the agent about new inputs and targets for the model. @@ -147,10 +149,10 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= # TODO: should be a check per model if len(self.table) > 2: - # if n_before_tell % self.train_every != n_after_tell % self.train_every: - self._update_models(train=train_models, a_priori_hypers=hypers) + if int(self.n_last_trained / self.train_every) != int(len(self.table) / self.train_every): + self._construct_models(train=train_models, a_priori_hypers=hypers) - def _update_models(self, train=True, skew_dims=None, a_priori_hypers=None): + def _construct_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 inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) @@ -173,7 +175,7 @@ def _update_models(self, train=True, skew_dims=None, a_priori_hypers=None): likelihood = gpytorch.likelihoods.GaussianLikelihood( noise_constraint=gpytorch.constraints.Interval( - torch.tensor(1e-2).square(), + torch.tensor(1e-4).square(), torch.tensor(1 / obj.min_snr).square(), ), # noise_prior=gpytorch.priors.torch_priors.LogNormalPrior(loc=loc, scale=scale), @@ -293,14 +295,12 @@ 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) @@ -619,7 +619,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): """ @@ -687,6 +687,8 @@ 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.""" diff --git a/bloptools/bayesian/dofs.py b/bloptools/bayesian/dofs.py index a68651f..5a14d75 100644 --- a/bloptools/bayesian/dofs.py +++ b/bloptools/bayesian/dofs.py @@ -9,7 +9,7 @@ from ophyd import Signal, SignalRO DEFAULT_BOUNDS = (-5.0, +5.0) -DOF_FIELDS = ["name", "description", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] +DOF_FIELDS = ["description", "readback", "lower_limit", "upper_limit", "units", "active", "read_only", "tags"] numeric = Union[float, int] @@ -110,9 +110,9 @@ def __repr__(self): @property def summary(self) -> pd.DataFrame: table = pd.DataFrame(columns=DOF_FIELDS) - for i, dof in enumerate(self.dofs): + for dof in self.dofs: for attr in table.columns: - table.loc[i, attr] = getattr(dof, attr) + table.loc[dof.name, attr] = getattr(dof, attr) # convert dtypes for attr in ["readback", "lower_limit", "upper_limit"]: diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index 7f15752..cc9842c 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -7,8 +7,8 @@ numeric = Union[float, int] -DEFAULT_MINIMUM_SNR = 1e1 -OBJ_FIELDS = ["description", "target", "limits", "weight", "log", "n", "snr", "min_snr"] +DEFAULT_MINIMUM_SNR = 1e2 +OBJ_FIELDS = ["description", "target", "active", "limits", "weight", "log", "n", "snr", "min_snr"] class DuplicateNameError(ValueError): @@ -30,6 +30,7 @@ class Objective: 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 diff --git a/bloptools/bayesian/plotting.py b/bloptools/bayesian/plotting.py index e774c95..e51256e 100644 --- a/bloptools/bayesian/plotting.py +++ b/bloptools/bayesian/plotting.py @@ -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"]) diff --git a/bloptools/tests/conftest.py b/bloptools/tests/conftest.py index 0fad97d..034b15c 100644 --- a/bloptools/tests/conftest.py +++ b/bloptools/tests/conftest.py @@ -71,7 +71,7 @@ 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): From 6364b6dd4fd4f002081e25b3491a6bdbbfa38ec4 Mon Sep 17 00:00:00 2001 From: Tanny Date: Mon, 20 Nov 2023 07:05:26 -0800 Subject: [PATCH 09/17] wrap up als work --- bloptools/bayesian/agent.py | 12 ++++++++++-- bloptools/bayesian/dofs.py | 5 ++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 8c6dbb9..e42517e 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -158,7 +158,9 @@ def _construct_models(self, train=True, skew_dims=None, a_priori_hypers=None): inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) for i, obj in enumerate(self.objectives): + values = self.get_objective_targets(i) + values = np.where(self.all_objectives_valid, values, np.nan) train_index = ~np.isnan(values) @@ -218,7 +220,7 @@ def _construct_models(self, train=True, skew_dims=None, a_priori_hypers=None): 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 ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True, upsample=1): """Ask the agent for the best point to sample, given an acquisition function. Parameters @@ -306,6 +308,11 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True): routing_index = utils.route(self.dofs.subset(active=True, read_only=False).readback, acq_points) acq_points = acq_points[routing_index] + 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) + return acq_points, acq_func_meta def acquire(self, acquisition_inputs): @@ -409,7 +416,7 @@ def learn( for i in range(iterations): print(f"running iteration {i + 1} / {iterations}") for single_acq_func in np.atleast_1d(acq_func): - acq_points, acq_func_meta = self.ask(n=n, acq_func_identifier=single_acq_func) + acq_points, acq_func_meta = self.ask(n=n, acq_func_identifier=single_acq_func, upsample=upsample) new_table = yield from self.acquire(acq_points) new_table.loc[:, "acq_func"] = acq_func_meta["name"] @@ -746,6 +753,7 @@ def go_to(self, **positions): def go_to_best(self): """Go to the position of the best input seen so far.""" yield from self.go_to(**self.best_inputs) + def plot_objectives(self, axes: Tuple = (0, 1), **kwargs): """Plot the sampled objectives diff --git a/bloptools/bayesian/dofs.py b/bloptools/bayesian/dofs.py index 5a14d75..e03764f 100644 --- a/bloptools/bayesian/dofs.py +++ b/bloptools/bayesian/dofs.py @@ -96,7 +96,10 @@ def __init__(self, dofs: list = []): self.dofs = dofs def __getitem__(self, i): - return self.dofs[i] + if type(i) is int: + return self.dofs[i] + elif type(i) is str: + return self.dofs[self.names.index(i)] def __len__(self): return len(self.dofs) From 82e98bdc36fc8a34359da13fdb709d3fd40ddbd8 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Mon, 20 Nov 2023 10:35:35 -0500 Subject: [PATCH 10/17] apply pre-commit --- bloptools/bayesian/agent.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index e42517e..611fdbf 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -158,7 +158,6 @@ def _construct_models(self, train=True, skew_dims=None, a_priori_hypers=None): inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) for i, obj in enumerate(self.objectives): - values = self.get_objective_targets(i) values = np.where(self.all_objectives_valid, values, np.nan) @@ -753,7 +752,6 @@ def go_to(self, **positions): def go_to_best(self): """Go to the position of the best input seen so far.""" yield from self.go_to(**self.best_inputs) - def plot_objectives(self, axes: Tuple = (0, 1), **kwargs): """Plot the sampled objectives From 51265863ffbb7ea3b8b5c4dbf178ea14c498fff8 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 28 Nov 2023 14:10:08 -0500 Subject: [PATCH 11/17] syntax for .tell() method is more flexible --- bloptools/bayesian/agent.py | 26 +- docs/source/tutorials/himmelblau.ipynb | 20 +- plot_benchmarks.ipynb | 354 ++++++++++++++++++++ run_benchmarks.ipynb | 432 +++++++++++++++++++++++++ 4 files changed, 827 insertions(+), 5 deletions(-) create mode 100644 plot_benchmarks.ipynb create mode 100644 run_benchmarks.ipynb diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 611fdbf..403caf4 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 @@ -119,7 +119,16 @@ def __init__( self.n_last_trained = 0 - def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models=True, hypers=None): + def tell( + self, + data: Optional[Mapping] = {}, + x: Optional[Mapping] = {}, + y: Optional[Mapping] = {}, + metadata: Optional[Mapping] = {}, + append=True, + train_models=True, + hypers=None, + ): """ Inform the agent about new inputs and targets for the model. @@ -139,9 +148,18 @@ def tell(self, x: Mapping, y: Mapping, metadata=None, append=True, train_models= A dict of hyperparameters for the model to assume a priori. """ - # n_before_tell = len(self.table) + 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({**x, **y, **metadata} if metadata is not None else {**x, **y}) + 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)) diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/himmelblau.ipynb index 29ceb16..802befa 100644 --- a/docs/source/tutorials/himmelblau.ipynb +++ b/docs/source/tutorials/himmelblau.ipynb @@ -170,6 +170,24 @@ "agent.plot_objectives()" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "75a6a54f-0140-4bd0-8367-0dd14cdf8116", + "metadata": {}, + "outputs": [], + "source": [ + "{x % 2 for x in np.arange(10)}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad31b676-8236-4cf8-b9be-c1214d2e8458", + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "id": "dc264346-10fb-4c88-9925-4bfcf0dd3b07", @@ -295,7 +313,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.10.0 ('nsls2')", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, diff --git a/plot_benchmarks.ipynb b/plot_benchmarks.ipynb new file mode 100644 index 0000000..b64ec04 --- /dev/null +++ b/plot_benchmarks.ipynb @@ -0,0 +1,354 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d778b648-f5b4-48a1-9c54-543f779f0271", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib as mpl\n", + "from matplotlib import pyplot as plt\n", + "from bloptools import test_functions\n", + "\n", + "import os\n", + "\n", + "os.environ[\n", + " \"PATH\"\n", + "] = \"/opt/homebrew/opt/llvm/bin:/Users/tom/opt/anaconda3/bin:/Users/tom/opt/anaconda3/condabin:/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Library/TeX/texbin:/opt/X11/bin\"\n", + "\n", + "plt.rcParams.update(\n", + " {\n", + " \"text.usetex\": True,\n", + " \"font.family\": \"serif\",\n", + " }\n", + ")\n", + "\n", + "cummax = lambda iterable: np.array([np.nanmax(iterable[:i]) for i in range(1, len(iterable) + 1)])\n", + "\n", + "\n", + "no_latent_paths = glob.glob(\"../data/benchmark_hartmann6_no_latent/*.h5\")\n", + "print(len(no_latent_paths))\n", + "no_latent_table = pd.DataFrame()\n", + "\n", + "for path in no_latent_paths:\n", + " no_latent_table.loc[:, path] = -cummax(np.exp(pd.read_hdf(path, key=\"table\").neg_hartmann_fitness.values))\n", + "\n", + "\n", + "latent_paths = sorted(glob.glob(\"../data/benchmark_hartmann6_latent-2/*.h5\"))\n", + "print(len(latent_paths))\n", + "latent_table = pd.DataFrame()\n", + "\n", + "for path in latent_paths:\n", + " # print(path)\n", + " latent_table.loc[:, path] = -cummax(np.exp(pd.read_hdf(path, key=\"table\").neg_hartmann_fitness.values))\n", + "\n", + "\n", + "pd.read_hdf(path, key=\"table\").acq_func.values[:10]\n", + "\n", + "max_iter = 256\n", + "\n", + "iter_bins = np.linspace(0, max_iter, 64)\n", + "\n", + "\n", + "import pandas as pd\n", + "import glob\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n", + "\n", + "\n", + "for col in no_latent_table.columns:\n", + " loss = no_latent_table.loc[:, col].values\n", + " ax.plot(np.arange(len(loss)) + 1, loss, lw=1e-1, c=\"r\")\n", + "\n", + "for col in latent_table.columns:\n", + " loss = latent_table.loc[:, col].values\n", + " ax.plot(np.arange(len(loss)) + 1, loss, lw=1e-1, c=\"b\")\n", + "\n", + "ax.plot(np.median(no_latent_table, axis=1), c=\"r\", lw=1e0)\n", + "ax.plot(np.median(latent_table, axis=1), c=\"b\", lw=1e0)\n", + "\n", + "xmin, xmax = 0, max_iter\n", + "ymin, ymax = ax.get_ylim()\n", + "\n", + "ax.set_xlabel(\"iteration\")\n", + "ax.set_ylabel(\"cumulative maximum\")\n", + "\n", + "optimum_inputs = None\n", + "optimum_target = -3.32237\n", + "\n", + "ax.plot([xmin, xmax], [optimum_target, optimum_target], ls=\":\", c=\"k\")\n", + "\n", + "ax.fill_betweenx([ymin, ymax], xmin, 64, alpha=0.1, label=\"quasi-random initialization\")\n", + "\n", + "ax.set_xlim(xmin, xmax)\n", + "ax.set_ylim(ymin, ymax)\n", + "\n", + "# ax.set_yscale(\"log\")\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6430af55-4737-4e3a-9652-d29a7632726f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a75629d0-0fa9-4153-8ea4-0d6caa29b8cc", + "metadata": {}, + "outputs": [], + "source": [ + "RE(agent.learn(\"qei\", n=2, iterations=16))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5767dbd1-6388-4da3-a33a-197142a6ba6e", + "metadata": {}, + "outputs": [], + "source": [ + "np.round(agent.tasks[0][\"model\"].covar_module.latent_transform.detach(), 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a42b5ec-1a55-4e23-9b0a-050e383d5349", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "405e21e2-a87d-4f79-8fb4-3bce7de37fce", + "metadata": {}, + "outputs": [], + "source": [ + "agent.plot_tasks()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f578fbad-03c1-4171-b27c-c266297443fc", + "metadata": {}, + "outputs": [], + "source": [ + "qei, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qei\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c582af1-02fa-4fc0-bfdb-16dfb525d5db", + "metadata": {}, + "outputs": [], + "source": [ + "qucb, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qucb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c32e76c-13db-4046-857c-6da522b94bc0", + "metadata": {}, + "outputs": [], + "source": [ + "X, _ = agent.ask(\"qei\", n=8)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb739668-9cb7-428d-bd63-a4cfa13fe0a1", + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "\n", + "x1 = torch.linspace(-6, 6, 63)\n", + "x2 = torch.linspace(-6, 6, 63)\n", + "X1, X2 = torch.meshgrid(x1, x2, indexing=\"ij\")\n", + "\n", + "xg = torch.cat([X1.unsqueeze(-1), X2.unsqueeze(-1)], dim=-1)\n", + "obj_grid = qei(xg.reshape(-1, 1, 2)).reshape(xg.shape[:2]).detach()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2e8f8cc-7f06-4fd8-b16b-e959b7fff934", + "metadata": {}, + "outputs": [], + "source": [ + "import scipy as sp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57f35ce5-3d4a-48d1-9ffe-750a35b26be8", + "metadata": {}, + "outputs": [], + "source": [ + "optima = [(3, 2), (-2.805, 3.1313), (-3.779, -3.283), (3.584, -1.848)]\n", + "\n", + "\n", + "fig, axes = plt.subplots(1, 4, figsize=(12, 4), dpi=256)\n", + "axes = axes.ravel()\n", + "\n", + "y = -agent._get_task_fitness(0)\n", + "\n", + "post = agent.model.posterior(xg)\n", + "pred = -post.mean[..., 0].detach()\n", + "err = post.variance[..., 0].detach()\n", + "\n", + "norm = mpl.colors.Normalize(-500, 0)\n", + "\n", + "cmap = \"magma\"\n", + "\n", + "zoom = 8\n", + "\n", + "uxg1 = sp.ndimage.zoom(xg[..., 0], zoom=zoom)\n", + "uxg2 = sp.ndimage.zoom(xg[..., 1], zoom=zoom)\n", + "umu = sp.ndimage.zoom(pred, zoom=zoom)\n", + "ulv = sp.ndimage.zoom(np.log(err.sqrt()), zoom=zoom)\n", + "uobj = sp.ndimage.zoom(obj_grid, zoom=zoom)\n", + "\n", + "samp = axes[0].scatter(*agent.active_inputs.values.T, c=-y, s=16, norm=norm, cmap=cmap)\n", + "# axes[1].pcolormesh(uxg1, uxg2, -umu, norm=norm, cmap=cmap)\n", + "# err = axes[2].pcolormesh(uxg1, uxg2, ulv, cmap=cmap)\n", + "# acq = axes[3].pcolormesh(uxg1, uxg2, uobj, cmap=cmap)\n", + "\n", + "\n", + "mean = axes[1].imshow(-umu.T[::-1], norm=norm, extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "lvar = axes[2].imshow(ulv.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "acqf = axes[3].imshow(uobj.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "\n", + "for ax in axes:\n", + " for x, y in optima:\n", + " ax.scatter(x, y, facecolor=\"none\", edgecolor=\"k\", marker=\"o\", lw=5e-1, s=16)\n", + "\n", + " ax.set_xlim(-6, 6)\n", + " ax.set_ylim(-6, 6)\n", + " ax.set_xticks(np.arange(-6, 7, 2))\n", + " ax.set_yticks(np.arange(-6, 7, 2))\n", + " ax.set_xlabel(\"$x_1$\")\n", + " ax.set_ylabel(\"$x_2$\")\n", + "\n", + "axes[3].scatter(*X.T, marker=\"x\", lw=5e-1, color=\"k\", s=16)\n", + "\n", + "# axes[0].set_title(\"sampled points\")\n", + "# axes[1].set_title(\"posterior mean\")\n", + "# axes[2].set_title(\"posterior log std. dev.\")\n", + "# axes[3].set_title(\"$q$-expected improvement\")\n", + "\n", + "\n", + "# clb1 = fig.colorbar(mappable=mpl.cm.ScalarMappable(norm=norm), ax=axes[0:2], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb1 = fig.colorbar(samp, ax=axes[0], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb1.set_label(\"sampled points\")\n", + "clb2 = fig.colorbar(mean, ax=axes[1], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb2.set_label(\"posterior mean\")\n", + "clb3 = fig.colorbar(lvar, ax=axes[2], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb3.set_label(\"posterior log std. dev.\")\n", + "clb4 = fig.colorbar(acqf, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb4.set_label(\"$q$-expected improvement\")\n", + "# clb3 = fig.colorbar(acq, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "966398d1-4d0b-4271-8b65-ac563e157347", + "metadata": {}, + "outputs": [], + "source": [ + "np.exp(ulv)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec8dacf5-2d2e-4f51-8c8b-83344d1465c5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0962910e-7c8f-4ccb-83ec-67551f6d443b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bda306b-6c3c-453e-8e1b-316c9001baee", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b897fa9e-22b7-4b0b-a57f-183cbf15bfc1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc343daa-643d-424f-a042-0e9ce7964adb", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(obj)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "543f65b5-7664-45ca-bd39-85d0b827febe", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/run_benchmarks.ipynb b/run_benchmarks.ipynb new file mode 100644 index 0000000..0144770 --- /dev/null +++ b/run_benchmarks.ipynb @@ -0,0 +1,432 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "7a05f2b4-ff97-4688-b9f6-e607ca4118c2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib as mpl\n", + "from matplotlib import pyplot as plt\n", + "from bloptools import test_functions\n", + "\n", + "import os\n", + "\n", + "os.environ[\n", + " \"PATH\"\n", + "] = \"/opt/homebrew/opt/llvm/bin:/Users/tom/opt/anaconda3/bin:/Users/tom/opt/anaconda3/condabin:/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Library/TeX/texbin:/opt/X11/bin\"\n", + "\n", + "plt.rcParams.update(\n", + " {\n", + " \"text.usetex\": True,\n", + " }\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef7272c7-9192-41a4-8b26-c620ddc2acfa", + "metadata": {}, + "outputs": [], + "source": [ + "from bloptools import devices\n", + "\n", + "dofs = [\n", + " {\"device\": devices.DOF(name=\"x1\"), \"limits\": (-6, 6), \"kind\": \"active\"},\n", + " {\"device\": devices.DOF(name=\"x2\"), \"limits\": (-6, 6), \"kind\": \"active\"},\n", + "]\n", + "\n", + "\n", + "def digestion(db, uid):\n", + " products = db[uid].table()\n", + "\n", + " for index, entry in products.iterrows():\n", + " products.loc[index, \"himmelblau\"] = test_functions.himmelblau(entry.x1, entry.x2)\n", + "\n", + " return products\n", + "\n", + "\n", + "tasks = [\n", + " {\"key\": \"himmelblau\", \"kind\": \"minimize\"},\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c39e262-6cc5-4165-bb06-16ce95641b14", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f315b273-8f5b-4cc8-93cd-317697e19841", + "metadata": {}, + "outputs": [], + "source": [ + "from bloptools import devices\n", + "\n", + "dofs = [\n", + " {\"device\": devices.DOF(name=\"x1\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + " {\"device\": devices.DOF(name=\"x2\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + " {\"device\": devices.DOF(name=\"x3\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + " {\"device\": devices.DOF(name=\"x4\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + " {\"device\": devices.DOF(name=\"x5\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + " {\"device\": devices.DOF(name=\"x6\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", + "]\n", + "\n", + "\n", + "def digestion(db, uid):\n", + " products = db[uid].table()\n", + "\n", + " for index, entry in products.iterrows():\n", + " products.loc[index, \"neg_hartmann\"] = -test_functions.hartmann6(\n", + " entry.x1, entry.x2, entry.x3, entry.x4, entry.x5, entry.x6\n", + " )\n", + "\n", + " return products\n", + "\n", + "\n", + "tasks = [\n", + " {\"key\": \"neg_hartmann\", \"kind\": \"maximize\", \"transform\": \"log\"},\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57e74305-9b50-4021-be84-07aea2ef1440", + "metadata": {}, + "outputs": [], + "source": [ + "import bloptools\n", + "from bloptools.utils import prepare_re_env\n", + "\n", + "%run -i $prepare_re_env.__file__ --db-type=temp\n", + "from bloptools.bayesian import Agent\n", + "\n", + "agent = Agent(\n", + " dofs=dofs,\n", + " tasks=tasks,\n", + " digestion=digestion,\n", + " db=db,\n", + ")\n", + "\n", + "# RE(agent.learn(\"qr\", n=128))\n", + "# RE(agent.learn(\"qei\", n=1, iterations=4))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae6d0836-7bed-4b69-a5d5-16fe2b00700a", + "metadata": {}, + "outputs": [], + "source": [ + "RE(\n", + " agent.benchmark(\n", + " output_dir=\"../data/benchmark_hartmann6_no_latent/\",\n", + " n_init=4,\n", + " runs=64,\n", + " learning_kwargs_list=[{\"acq_func\": \"qei\", \"n\": 4, \"iterations\": 64}],\n", + " )\n", + ")\n", + "# {\"acq_func\": \"qem\", \"n\": 4, \"iterations\": 16},\n", + "# {\"acq_func\": \"qei\", \"n\": 4, \"iterations\": 16},\n", + "# {\"acq_func\": \"qem\", \"n\": 4, \"iterations\": 16}]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e35d047-7354-478b-8bb6-55d39f7727d9", + "metadata": {}, + "outputs": [], + "source": [ + "os.mkdir()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97d534dd-af7d-46db-af7c-2464a2da119e", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import glob\n", + "\n", + "paths = glob.glob(\"../data/benchmark_hartmann_no_latent/*.h5\")\n", + "\n", + "for path in paths:\n", + " table = pd.read_hdf(path, key=\"table\")\n", + "\n", + " plt.plot(-cummax(np.exp(table.neg_hartmann_fitness.values)), lw=1e-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6430af55-4737-4e3a-9652-d29a7632726f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a75629d0-0fa9-4153-8ea4-0d6caa29b8cc", + "metadata": {}, + "outputs": [], + "source": [ + "RE(agent.learn(\"qei\", n=2, iterations=16))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5767dbd1-6388-4da3-a33a-197142a6ba6e", + "metadata": {}, + "outputs": [], + "source": [ + "np.round(agent.tasks[0][\"model\"].covar_module.latent_transform.detach(), 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a42b5ec-1a55-4e23-9b0a-050e383d5349", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7f1e3cc-0d93-4c35-88e9-47811f2a5fb7", + "metadata": {}, + "outputs": [], + "source": [ + "cummax = lambda iterable: np.array([np.nanmax(iterable[:i]) for i in range(1, len(iterable) + 1)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "405e21e2-a87d-4f79-8fb4-3bce7de37fce", + "metadata": {}, + "outputs": [], + "source": [ + "agent.plot_tasks()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f578fbad-03c1-4171-b27c-c266297443fc", + "metadata": {}, + "outputs": [], + "source": [ + "qei, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qei\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c32e76c-13db-4046-857c-6da522b94bc0", + "metadata": {}, + "outputs": [], + "source": [ + "X, _ = agent.ask(\"qei\", n=8)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb739668-9cb7-428d-bd63-a4cfa13fe0a1", + "metadata": {}, + "outputs": [], + "source": [ + "import torch\n", + "\n", + "x1 = torch.linspace(-6, 6, 63)\n", + "x2 = torch.linspace(-6, 6, 63)\n", + "X1, X2 = torch.meshgrid(x1, x2, indexing=\"ij\")\n", + "\n", + "xg = torch.cat([X1.unsqueeze(-1), X2.unsqueeze(-1)], dim=-1)\n", + "obj_grid = qei(xg.reshape(-1, 1, 2)).reshape(xg.shape[:2]).detach()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2e8f8cc-7f06-4fd8-b16b-e959b7fff934", + "metadata": {}, + "outputs": [], + "source": [ + "import scipy as sp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57f35ce5-3d4a-48d1-9ffe-750a35b26be8", + "metadata": {}, + "outputs": [], + "source": [ + "optima = [(3, 2), (-2.805, 3.1313), (-3.779, -3.283), (3.584, -1.848)]\n", + "\n", + "\n", + "fig, axes = plt.subplots(1, 4, figsize=(12, 4), dpi=256)\n", + "axes = axes.ravel()\n", + "\n", + "y = -agent._get_task_fitness(0)\n", + "\n", + "post = agent.model.posterior(xg)\n", + "pred = -post.mean[..., 0].detach()\n", + "err = post.variance[..., 0].detach()\n", + "\n", + "norm = mpl.colors.Normalize(-500, 0)\n", + "\n", + "cmap = \"magma\"\n", + "\n", + "zoom = 8\n", + "\n", + "uxg1 = sp.ndimage.zoom(xg[..., 0], zoom=zoom)\n", + "uxg2 = sp.ndimage.zoom(xg[..., 1], zoom=zoom)\n", + "umu = sp.ndimage.zoom(pred, zoom=zoom)\n", + "usig = np.exp(sp.ndimage.zoom(np.log(err.sqrt()), zoom=zoom))\n", + "uobj = sp.ndimage.zoom(obj_grid, zoom=zoom)\n", + "\n", + "samp = axes[0].scatter(*agent.active_inputs.values.T, c=-y, s=16, norm=norm, cmap=cmap)\n", + "# axes[1].pcolormesh(uxg1, uxg2, -umu, norm=norm, cmap=cmap)\n", + "# err = axes[2].pcolormesh(uxg1, uxg2, ulv, cmap=cmap)\n", + "# acq = axes[3].pcolormesh(uxg1, uxg2, uobj, cmap=cmap)\n", + "\n", + "\n", + "mean = axes[1].imshow(-umu.T[::-1], norm=norm, extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "err = axes[2].imshow(usig.T[::-1], norm=mpl.colors.LogNorm(), extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "acqf = axes[3].imshow(uobj.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", + "\n", + "for ax in axes:\n", + " for x, y in optima:\n", + " ax.scatter(x, y, facecolor=\"w\", edgecolor=\"k\", marker=\"o\", lw=5e-1, s=16)\n", + "\n", + " ax.set_xlim(-6, 6)\n", + " ax.set_ylim(-6, 6)\n", + " ax.set_xticks(np.arange(-6, 7, 2))\n", + " ax.set_yticks(np.arange(-6, 7, 2))\n", + " ax.set_xlabel(\"$x_1$\")\n", + " ax.set_ylabel(\"$x_2$\")\n", + "\n", + "axes[3].scatter(*X.T, marker=\"d\", lw=5e-1, color=\"k\", s=16)\n", + "\n", + "# axes[0].set_title(\"sampled points\")\n", + "# axes[1].set_title(\"posterior mean\")\n", + "# axes[2].set_title(\"posterior log std. dev.\")\n", + "# axes[3].set_title(\"$q$-expected improvement\")\n", + "\n", + "\n", + "# clb1 = fig.colorbar(mappable=mpl.cm.ScalarMappable(norm=norm), ax=axes[0:2], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb1 = fig.colorbar(samp, ax=axes[0], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb1.set_label(\"sampled points\")\n", + "clb2 = fig.colorbar(mean, ax=axes[1], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb2.set_label(\"posterior mean\")\n", + "clb3 = fig.colorbar(err, ax=axes[2], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb3.set_label(\"posterior std. dev.\")\n", + "clb4 = fig.colorbar(acqf, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", + "clb4.set_label(\"$q$-expected improvement\")\n", + "# clb3 = fig.colorbar(acq, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "966398d1-4d0b-4271-8b65-ac563e157347", + "metadata": {}, + "outputs": [], + "source": [ + "np.exp(ulv)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec8dacf5-2d2e-4f51-8c8b-83344d1465c5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0962910e-7c8f-4ccb-83ec-67551f6d443b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bda306b-6c3c-453e-8e1b-316c9001baee", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b897fa9e-22b7-4b0b-a57f-183cbf15bfc1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc343daa-643d-424f-a042-0e9ce7964adb", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(obj)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "543f65b5-7664-45ca-bd39-85d0b827febe", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 756868366b7a4abf34c60153d287007f2a7bba35 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 28 Nov 2023 14:10:42 -0500 Subject: [PATCH 12/17] remove stray notebooks --- plot_benchmarks.ipynb | 354 ---------------------------------- run_benchmarks.ipynb | 432 ------------------------------------------ 2 files changed, 786 deletions(-) delete mode 100644 plot_benchmarks.ipynb delete mode 100644 run_benchmarks.ipynb diff --git a/plot_benchmarks.ipynb b/plot_benchmarks.ipynb deleted file mode 100644 index b64ec04..0000000 --- a/plot_benchmarks.ipynb +++ /dev/null @@ -1,354 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "d778b648-f5b4-48a1-9c54-543f779f0271", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib as mpl\n", - "from matplotlib import pyplot as plt\n", - "from bloptools import test_functions\n", - "\n", - "import os\n", - "\n", - "os.environ[\n", - " \"PATH\"\n", - "] = \"/opt/homebrew/opt/llvm/bin:/Users/tom/opt/anaconda3/bin:/Users/tom/opt/anaconda3/condabin:/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Library/TeX/texbin:/opt/X11/bin\"\n", - "\n", - "plt.rcParams.update(\n", - " {\n", - " \"text.usetex\": True,\n", - " \"font.family\": \"serif\",\n", - " }\n", - ")\n", - "\n", - "cummax = lambda iterable: np.array([np.nanmax(iterable[:i]) for i in range(1, len(iterable) + 1)])\n", - "\n", - "\n", - "no_latent_paths = glob.glob(\"../data/benchmark_hartmann6_no_latent/*.h5\")\n", - "print(len(no_latent_paths))\n", - "no_latent_table = pd.DataFrame()\n", - "\n", - "for path in no_latent_paths:\n", - " no_latent_table.loc[:, path] = -cummax(np.exp(pd.read_hdf(path, key=\"table\").neg_hartmann_fitness.values))\n", - "\n", - "\n", - "latent_paths = sorted(glob.glob(\"../data/benchmark_hartmann6_latent-2/*.h5\"))\n", - "print(len(latent_paths))\n", - "latent_table = pd.DataFrame()\n", - "\n", - "for path in latent_paths:\n", - " # print(path)\n", - " latent_table.loc[:, path] = -cummax(np.exp(pd.read_hdf(path, key=\"table\").neg_hartmann_fitness.values))\n", - "\n", - "\n", - "pd.read_hdf(path, key=\"table\").acq_func.values[:10]\n", - "\n", - "max_iter = 256\n", - "\n", - "iter_bins = np.linspace(0, max_iter, 64)\n", - "\n", - "\n", - "import pandas as pd\n", - "import glob\n", - "\n", - "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n", - "\n", - "\n", - "for col in no_latent_table.columns:\n", - " loss = no_latent_table.loc[:, col].values\n", - " ax.plot(np.arange(len(loss)) + 1, loss, lw=1e-1, c=\"r\")\n", - "\n", - "for col in latent_table.columns:\n", - " loss = latent_table.loc[:, col].values\n", - " ax.plot(np.arange(len(loss)) + 1, loss, lw=1e-1, c=\"b\")\n", - "\n", - "ax.plot(np.median(no_latent_table, axis=1), c=\"r\", lw=1e0)\n", - "ax.plot(np.median(latent_table, axis=1), c=\"b\", lw=1e0)\n", - "\n", - "xmin, xmax = 0, max_iter\n", - "ymin, ymax = ax.get_ylim()\n", - "\n", - "ax.set_xlabel(\"iteration\")\n", - "ax.set_ylabel(\"cumulative maximum\")\n", - "\n", - "optimum_inputs = None\n", - "optimum_target = -3.32237\n", - "\n", - "ax.plot([xmin, xmax], [optimum_target, optimum_target], ls=\":\", c=\"k\")\n", - "\n", - "ax.fill_betweenx([ymin, ymax], xmin, 64, alpha=0.1, label=\"quasi-random initialization\")\n", - "\n", - "ax.set_xlim(xmin, xmax)\n", - "ax.set_ylim(ymin, ymax)\n", - "\n", - "# ax.set_yscale(\"log\")\n", - "\n", - "plt.tight_layout()\n", - "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6430af55-4737-4e3a-9652-d29a7632726f", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a75629d0-0fa9-4153-8ea4-0d6caa29b8cc", - "metadata": {}, - "outputs": [], - "source": [ - "RE(agent.learn(\"qei\", n=2, iterations=16))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5767dbd1-6388-4da3-a33a-197142a6ba6e", - "metadata": {}, - "outputs": [], - "source": [ - "np.round(agent.tasks[0][\"model\"].covar_module.latent_transform.detach(), 3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7a42b5ec-1a55-4e23-9b0a-050e383d5349", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "405e21e2-a87d-4f79-8fb4-3bce7de37fce", - "metadata": {}, - "outputs": [], - "source": [ - "agent.plot_tasks()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f578fbad-03c1-4171-b27c-c266297443fc", - "metadata": {}, - "outputs": [], - "source": [ - "qei, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qei\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5c582af1-02fa-4fc0-bfdb-16dfb525d5db", - "metadata": {}, - "outputs": [], - "source": [ - "qucb, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qucb\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4c32e76c-13db-4046-857c-6da522b94bc0", - "metadata": {}, - "outputs": [], - "source": [ - "X, _ = agent.ask(\"qei\", n=8)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fb739668-9cb7-428d-bd63-a4cfa13fe0a1", - "metadata": {}, - "outputs": [], - "source": [ - "import torch\n", - "\n", - "x1 = torch.linspace(-6, 6, 63)\n", - "x2 = torch.linspace(-6, 6, 63)\n", - "X1, X2 = torch.meshgrid(x1, x2, indexing=\"ij\")\n", - "\n", - "xg = torch.cat([X1.unsqueeze(-1), X2.unsqueeze(-1)], dim=-1)\n", - "obj_grid = qei(xg.reshape(-1, 1, 2)).reshape(xg.shape[:2]).detach()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e2e8f8cc-7f06-4fd8-b16b-e959b7fff934", - "metadata": {}, - "outputs": [], - "source": [ - "import scipy as sp" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "57f35ce5-3d4a-48d1-9ffe-750a35b26be8", - "metadata": {}, - "outputs": [], - "source": [ - "optima = [(3, 2), (-2.805, 3.1313), (-3.779, -3.283), (3.584, -1.848)]\n", - "\n", - "\n", - "fig, axes = plt.subplots(1, 4, figsize=(12, 4), dpi=256)\n", - "axes = axes.ravel()\n", - "\n", - "y = -agent._get_task_fitness(0)\n", - "\n", - "post = agent.model.posterior(xg)\n", - "pred = -post.mean[..., 0].detach()\n", - "err = post.variance[..., 0].detach()\n", - "\n", - "norm = mpl.colors.Normalize(-500, 0)\n", - "\n", - "cmap = \"magma\"\n", - "\n", - "zoom = 8\n", - "\n", - "uxg1 = sp.ndimage.zoom(xg[..., 0], zoom=zoom)\n", - "uxg2 = sp.ndimage.zoom(xg[..., 1], zoom=zoom)\n", - "umu = sp.ndimage.zoom(pred, zoom=zoom)\n", - "ulv = sp.ndimage.zoom(np.log(err.sqrt()), zoom=zoom)\n", - "uobj = sp.ndimage.zoom(obj_grid, zoom=zoom)\n", - "\n", - "samp = axes[0].scatter(*agent.active_inputs.values.T, c=-y, s=16, norm=norm, cmap=cmap)\n", - "# axes[1].pcolormesh(uxg1, uxg2, -umu, norm=norm, cmap=cmap)\n", - "# err = axes[2].pcolormesh(uxg1, uxg2, ulv, cmap=cmap)\n", - "# acq = axes[3].pcolormesh(uxg1, uxg2, uobj, cmap=cmap)\n", - "\n", - "\n", - "mean = axes[1].imshow(-umu.T[::-1], norm=norm, extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "lvar = axes[2].imshow(ulv.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "acqf = axes[3].imshow(uobj.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "\n", - "for ax in axes:\n", - " for x, y in optima:\n", - " ax.scatter(x, y, facecolor=\"none\", edgecolor=\"k\", marker=\"o\", lw=5e-1, s=16)\n", - "\n", - " ax.set_xlim(-6, 6)\n", - " ax.set_ylim(-6, 6)\n", - " ax.set_xticks(np.arange(-6, 7, 2))\n", - " ax.set_yticks(np.arange(-6, 7, 2))\n", - " ax.set_xlabel(\"$x_1$\")\n", - " ax.set_ylabel(\"$x_2$\")\n", - "\n", - "axes[3].scatter(*X.T, marker=\"x\", lw=5e-1, color=\"k\", s=16)\n", - "\n", - "# axes[0].set_title(\"sampled points\")\n", - "# axes[1].set_title(\"posterior mean\")\n", - "# axes[2].set_title(\"posterior log std. dev.\")\n", - "# axes[3].set_title(\"$q$-expected improvement\")\n", - "\n", - "\n", - "# clb1 = fig.colorbar(mappable=mpl.cm.ScalarMappable(norm=norm), ax=axes[0:2], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb1 = fig.colorbar(samp, ax=axes[0], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb1.set_label(\"sampled points\")\n", - "clb2 = fig.colorbar(mean, ax=axes[1], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb2.set_label(\"posterior mean\")\n", - "clb3 = fig.colorbar(lvar, ax=axes[2], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb3.set_label(\"posterior log std. dev.\")\n", - "clb4 = fig.colorbar(acqf, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb4.set_label(\"$q$-expected improvement\")\n", - "# clb3 = fig.colorbar(acq, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", - "\n", - "plt.tight_layout()\n", - "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "966398d1-4d0b-4271-8b65-ac563e157347", - "metadata": {}, - "outputs": [], - "source": [ - "np.exp(ulv)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ec8dacf5-2d2e-4f51-8c8b-83344d1465c5", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0962910e-7c8f-4ccb-83ec-67551f6d443b", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5bda306b-6c3c-453e-8e1b-316c9001baee", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b897fa9e-22b7-4b0b-a57f-183cbf15bfc1", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cc343daa-643d-424f-a042-0e9ce7964adb", - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(obj)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "543f65b5-7664-45ca-bd39-85d0b827febe", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.0" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/run_benchmarks.ipynb b/run_benchmarks.ipynb deleted file mode 100644 index 0144770..0000000 --- a/run_benchmarks.ipynb +++ /dev/null @@ -1,432 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "7a05f2b4-ff97-4688-b9f6-e607ca4118c2", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib as mpl\n", - "from matplotlib import pyplot as plt\n", - "from bloptools import test_functions\n", - "\n", - "import os\n", - "\n", - "os.environ[\n", - " \"PATH\"\n", - "] = \"/opt/homebrew/opt/llvm/bin:/Users/tom/opt/anaconda3/bin:/Users/tom/opt/anaconda3/condabin:/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin:/Library/TeX/texbin:/opt/X11/bin\"\n", - "\n", - "plt.rcParams.update(\n", - " {\n", - " \"text.usetex\": True,\n", - " }\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ef7272c7-9192-41a4-8b26-c620ddc2acfa", - "metadata": {}, - "outputs": [], - "source": [ - "from bloptools import devices\n", - "\n", - "dofs = [\n", - " {\"device\": devices.DOF(name=\"x1\"), \"limits\": (-6, 6), \"kind\": \"active\"},\n", - " {\"device\": devices.DOF(name=\"x2\"), \"limits\": (-6, 6), \"kind\": \"active\"},\n", - "]\n", - "\n", - "\n", - "def digestion(db, uid):\n", - " products = db[uid].table()\n", - "\n", - " for index, entry in products.iterrows():\n", - " products.loc[index, \"himmelblau\"] = test_functions.himmelblau(entry.x1, entry.x2)\n", - "\n", - " return products\n", - "\n", - "\n", - "tasks = [\n", - " {\"key\": \"himmelblau\", \"kind\": \"minimize\"},\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5c39e262-6cc5-4165-bb06-16ce95641b14", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f315b273-8f5b-4cc8-93cd-317697e19841", - "metadata": {}, - "outputs": [], - "source": [ - "from bloptools import devices\n", - "\n", - "dofs = [\n", - " {\"device\": devices.DOF(name=\"x1\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - " {\"device\": devices.DOF(name=\"x2\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - " {\"device\": devices.DOF(name=\"x3\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - " {\"device\": devices.DOF(name=\"x4\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - " {\"device\": devices.DOF(name=\"x5\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - " {\"device\": devices.DOF(name=\"x6\"), \"limits\": (0, 1), \"kind\": \"active\"}, # \"latent_group\": 0},\n", - "]\n", - "\n", - "\n", - "def digestion(db, uid):\n", - " products = db[uid].table()\n", - "\n", - " for index, entry in products.iterrows():\n", - " products.loc[index, \"neg_hartmann\"] = -test_functions.hartmann6(\n", - " entry.x1, entry.x2, entry.x3, entry.x4, entry.x5, entry.x6\n", - " )\n", - "\n", - " return products\n", - "\n", - "\n", - "tasks = [\n", - " {\"key\": \"neg_hartmann\", \"kind\": \"maximize\", \"transform\": \"log\"},\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "57e74305-9b50-4021-be84-07aea2ef1440", - "metadata": {}, - "outputs": [], - "source": [ - "import bloptools\n", - "from bloptools.utils import prepare_re_env\n", - "\n", - "%run -i $prepare_re_env.__file__ --db-type=temp\n", - "from bloptools.bayesian import Agent\n", - "\n", - "agent = Agent(\n", - " dofs=dofs,\n", - " tasks=tasks,\n", - " digestion=digestion,\n", - " db=db,\n", - ")\n", - "\n", - "# RE(agent.learn(\"qr\", n=128))\n", - "# RE(agent.learn(\"qei\", n=1, iterations=4))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ae6d0836-7bed-4b69-a5d5-16fe2b00700a", - "metadata": {}, - "outputs": [], - "source": [ - "RE(\n", - " agent.benchmark(\n", - " output_dir=\"../data/benchmark_hartmann6_no_latent/\",\n", - " n_init=4,\n", - " runs=64,\n", - " learning_kwargs_list=[{\"acq_func\": \"qei\", \"n\": 4, \"iterations\": 64}],\n", - " )\n", - ")\n", - "# {\"acq_func\": \"qem\", \"n\": 4, \"iterations\": 16},\n", - "# {\"acq_func\": \"qei\", \"n\": 4, \"iterations\": 16},\n", - "# {\"acq_func\": \"qem\", \"n\": 4, \"iterations\": 16}]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2e35d047-7354-478b-8bb6-55d39f7727d9", - "metadata": {}, - "outputs": [], - "source": [ - "os.mkdir()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "97d534dd-af7d-46db-af7c-2464a2da119e", - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "import glob\n", - "\n", - "paths = glob.glob(\"../data/benchmark_hartmann_no_latent/*.h5\")\n", - "\n", - "for path in paths:\n", - " table = pd.read_hdf(path, key=\"table\")\n", - "\n", - " plt.plot(-cummax(np.exp(table.neg_hartmann_fitness.values)), lw=1e-1)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6430af55-4737-4e3a-9652-d29a7632726f", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a75629d0-0fa9-4153-8ea4-0d6caa29b8cc", - "metadata": {}, - "outputs": [], - "source": [ - "RE(agent.learn(\"qei\", n=2, iterations=16))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5767dbd1-6388-4da3-a33a-197142a6ba6e", - "metadata": {}, - "outputs": [], - "source": [ - "np.round(agent.tasks[0][\"model\"].covar_module.latent_transform.detach(), 3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7a42b5ec-1a55-4e23-9b0a-050e383d5349", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e7f1e3cc-0d93-4c35-88e9-47811f2a5fb7", - "metadata": {}, - "outputs": [], - "source": [ - "cummax = lambda iterable: np.array([np.nanmax(iterable[:i]) for i in range(1, len(iterable) + 1)])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "405e21e2-a87d-4f79-8fb4-3bce7de37fce", - "metadata": {}, - "outputs": [], - "source": [ - "agent.plot_tasks()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f578fbad-03c1-4171-b27c-c266297443fc", - "metadata": {}, - "outputs": [], - "source": [ - "qei, _ = bloptools.bayesian.acquisition.get_acquisition_function(agent, \"qei\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4c32e76c-13db-4046-857c-6da522b94bc0", - "metadata": {}, - "outputs": [], - "source": [ - "X, _ = agent.ask(\"qei\", n=8)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fb739668-9cb7-428d-bd63-a4cfa13fe0a1", - "metadata": {}, - "outputs": [], - "source": [ - "import torch\n", - "\n", - "x1 = torch.linspace(-6, 6, 63)\n", - "x2 = torch.linspace(-6, 6, 63)\n", - "X1, X2 = torch.meshgrid(x1, x2, indexing=\"ij\")\n", - "\n", - "xg = torch.cat([X1.unsqueeze(-1), X2.unsqueeze(-1)], dim=-1)\n", - "obj_grid = qei(xg.reshape(-1, 1, 2)).reshape(xg.shape[:2]).detach()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e2e8f8cc-7f06-4fd8-b16b-e959b7fff934", - "metadata": {}, - "outputs": [], - "source": [ - "import scipy as sp" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "57f35ce5-3d4a-48d1-9ffe-750a35b26be8", - "metadata": {}, - "outputs": [], - "source": [ - "optima = [(3, 2), (-2.805, 3.1313), (-3.779, -3.283), (3.584, -1.848)]\n", - "\n", - "\n", - "fig, axes = plt.subplots(1, 4, figsize=(12, 4), dpi=256)\n", - "axes = axes.ravel()\n", - "\n", - "y = -agent._get_task_fitness(0)\n", - "\n", - "post = agent.model.posterior(xg)\n", - "pred = -post.mean[..., 0].detach()\n", - "err = post.variance[..., 0].detach()\n", - "\n", - "norm = mpl.colors.Normalize(-500, 0)\n", - "\n", - "cmap = \"magma\"\n", - "\n", - "zoom = 8\n", - "\n", - "uxg1 = sp.ndimage.zoom(xg[..., 0], zoom=zoom)\n", - "uxg2 = sp.ndimage.zoom(xg[..., 1], zoom=zoom)\n", - "umu = sp.ndimage.zoom(pred, zoom=zoom)\n", - "usig = np.exp(sp.ndimage.zoom(np.log(err.sqrt()), zoom=zoom))\n", - "uobj = sp.ndimage.zoom(obj_grid, zoom=zoom)\n", - "\n", - "samp = axes[0].scatter(*agent.active_inputs.values.T, c=-y, s=16, norm=norm, cmap=cmap)\n", - "# axes[1].pcolormesh(uxg1, uxg2, -umu, norm=norm, cmap=cmap)\n", - "# err = axes[2].pcolormesh(uxg1, uxg2, ulv, cmap=cmap)\n", - "# acq = axes[3].pcolormesh(uxg1, uxg2, uobj, cmap=cmap)\n", - "\n", - "\n", - "mean = axes[1].imshow(-umu.T[::-1], norm=norm, extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "err = axes[2].imshow(usig.T[::-1], norm=mpl.colors.LogNorm(), extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "acqf = axes[3].imshow(uobj.T[::-1], extent=[-6, 6, -6, 6], cmap=cmap, aspect=\"auto\")\n", - "\n", - "for ax in axes:\n", - " for x, y in optima:\n", - " ax.scatter(x, y, facecolor=\"w\", edgecolor=\"k\", marker=\"o\", lw=5e-1, s=16)\n", - "\n", - " ax.set_xlim(-6, 6)\n", - " ax.set_ylim(-6, 6)\n", - " ax.set_xticks(np.arange(-6, 7, 2))\n", - " ax.set_yticks(np.arange(-6, 7, 2))\n", - " ax.set_xlabel(\"$x_1$\")\n", - " ax.set_ylabel(\"$x_2$\")\n", - "\n", - "axes[3].scatter(*X.T, marker=\"d\", lw=5e-1, color=\"k\", s=16)\n", - "\n", - "# axes[0].set_title(\"sampled points\")\n", - "# axes[1].set_title(\"posterior mean\")\n", - "# axes[2].set_title(\"posterior log std. dev.\")\n", - "# axes[3].set_title(\"$q$-expected improvement\")\n", - "\n", - "\n", - "# clb1 = fig.colorbar(mappable=mpl.cm.ScalarMappable(norm=norm), ax=axes[0:2], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb1 = fig.colorbar(samp, ax=axes[0], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb1.set_label(\"sampled points\")\n", - "clb2 = fig.colorbar(mean, ax=axes[1], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb2.set_label(\"posterior mean\")\n", - "clb3 = fig.colorbar(err, ax=axes[2], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb3.set_label(\"posterior std. dev.\")\n", - "clb4 = fig.colorbar(acqf, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", - "clb4.set_label(\"$q$-expected improvement\")\n", - "# clb3 = fig.colorbar(acq, ax=axes[3], location=\"bottom\", shrink=0.8, aspect=32)\n", - "\n", - "plt.tight_layout()\n", - "plt.savefig(\"../plots/bayesian_himmelblau.pdf\", bbox_inches=\"tight\", transparent=True, pad_inches=0, dpi=256)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "966398d1-4d0b-4271-8b65-ac563e157347", - "metadata": {}, - "outputs": [], - "source": [ - "np.exp(ulv)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ec8dacf5-2d2e-4f51-8c8b-83344d1465c5", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0962910e-7c8f-4ccb-83ec-67551f6d443b", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5bda306b-6c3c-453e-8e1b-316c9001baee", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b897fa9e-22b7-4b0b-a57f-183cbf15bfc1", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cc343daa-643d-424f-a042-0e9ce7964adb", - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(obj)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "543f65b5-7664-45ca-bd39-85d0b827febe", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.0" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From 9aae8cc775be347fd0591ed321835e4b9e9a1254 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 28 Nov 2023 14:22:56 -0500 Subject: [PATCH 13/17] fix .reset() method --- bloptools/bayesian/agent.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 403caf4..8c9c08c 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -163,13 +163,15 @@ def tell( self.table = pd.concat([self.table, new_table]) if append else new_table self.table.index = np.arange(len(self.table)) - # n_after_tell = len(self.table) - # TODO: should be a check per model if len(self.table) > 2: - if int(self.n_last_trained / self.train_every) != int(len(self.table) / self.train_every): + if any([not hasattr(obj, "model") for obj in self.objectives]): self._construct_models(train=train_models, a_priori_hypers=hypers) + elif int(self.n_last_trained / self.train_every) != int(len(self.table) / self.train_every): + if train_models: + self._train_models(a_priori_hypers=hypers) + def _construct_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 @@ -451,8 +453,12 @@ def get_acquisition_function(self, identifier, return_metadata=False): def reset(self): """Reset the agent.""" self.table = pd.DataFrame() + for obj in self.objectives: - del obj.model + 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}] From f14bd5570e20ab9307ece14a46217959bedf7993 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Sat, 6 Jan 2024 15:57:31 -0500 Subject: [PATCH 14/17] addressed PR review --- bloptools/bayesian/__init__.py | 1 - bloptools/bayesian/devices.py | 59 -------------------------------- bloptools/bayesian/dofs.py | 50 ++++++++++++++++++--------- bloptools/bayesian/objectives.py | 42 ++++++++++++++++++----- 4 files changed, 68 insertions(+), 84 deletions(-) delete mode 100644 bloptools/bayesian/devices.py diff --git a/bloptools/bayesian/__init__.py b/bloptools/bayesian/__init__.py index 802690e..e49a76f 100644 --- a/bloptools/bayesian/__init__.py +++ b/bloptools/bayesian/__init__.py @@ -1,4 +1,3 @@ from .agent import * # noqa F401 -from .devices import * # noqa F401 from .dofs import * # noqa F401 from .objectives import * # noqa F401 diff --git a/bloptools/bayesian/devices.py b/bloptools/bayesian/devices.py deleted file mode 100644 index 3208a4d..0000000 --- a/bloptools/bayesian/devices.py +++ /dev/null @@ -1,59 +0,0 @@ -import time as ttime -import uuid - -import numpy as np -from ophyd import Signal, SignalRO # noqa F401 - - -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. - """ - - def __init__(self, name=None, theta=0.95, *args, **kwargs): - name = name if name is not None else str(uuid.uuid4()) - - super().__init__(name=name, *args, **kwargs) - - self.theta = theta - self.old_t = ttime.monotonic() - self.old_y = 0.0 - - def get(self): - new_t = ttime.monotonic() - alpha = self.theta ** (new_t - self.old_t) - new_y = alpha * self.old_y + np.sqrt(1 - alpha**2) * np.random.standard_normal() - - self.old_t = new_t - self.old_y = new_y - return new_y - - -class TimeReadback(SignalRO): - """ - Returns the current timestamp. - """ - - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - - def get(self): - return ttime.time() - - -class ConstantReadback(SignalRO): - """ - Returns a constant every time you read it (more useful than you'd think). - """ - - def __init__(self, constant=1, *args, **kwargs): - super().__init__(*args, **kwargs) - - self.constant = constant - - def get(self): - return self.constant diff --git a/bloptools/bayesian/dofs.py b/bloptools/bayesian/dofs.py index e03764f..ad3f39f 100644 --- a/bloptools/bayesian/dofs.py +++ b/bloptools/bayesian/dofs.py @@ -25,13 +25,42 @@ def _validate_dofs(dofs): 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. + + Parameters + ---------- + name: str + The name of the DOF. This is used as a key. + description: str + A longer name 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. + """ + device: Signal = None description: str = None name: str = None @@ -40,8 +69,10 @@ class DOF: 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()) @@ -100,6 +131,8 @@ def __getitem__(self, i): 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): return len(self.dofs) @@ -179,21 +212,6 @@ def _dof_mask(self, active=None, read_only=None, tags=[]): 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 _subset_devices(self, read_only=None, active=None, tags=[]): - # return [dof["device"] for dof in self._subset_dofs(read_only, active, tags)] - - # def _read_subset_devices(self, read_only=None, active=None, tags=[]): - # return [device.read()[device.name]["value"] for device in self._subset_devices(read_only, active, tags)] - - # def _subset_dof_names(self, read_only=None, active=None, tags=[]): - # return [device.name for device in self._subset_devices(read_only, active, tags)] - - # def _subset_dof_limits(self, read_only=None, active=None, tags=[]): - # dofs_subset = self._subset_dofs(read_only, active, tags) - # if len(dofs_subset) > 0: - # return torch.tensor([dof["limits"] for dof in dofs_subset], dtype=torch.float64).T - # return torch.empty((2, 0)) - def activate(self, read_only=None, active=None, tags=[]): for dof in self._subset_dofs(read_only, active, tags): dof.active = True diff --git a/bloptools/bayesian/objectives.py b/bloptools/bayesian/objectives.py index cc9842c..7bec712 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -20,11 +20,36 @@ def _validate_objectives(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}"') + 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" @@ -46,8 +71,6 @@ def __post_init__(self): if self.target not in ["min", "max"]: raise ValueError("'target' must be either 'min', 'max', or a number.") - # self.device = Signal(name=self.name) - @property def label(self): return f"{'log ' if self.log else ''}{self.description}" @@ -102,21 +125,24 @@ def summary(self): def __repr__(self): return self.summary.__repr__() - # @property - # def descriptions(self) -> list: - # return [obj.description for obj in self.objectives] + @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 weights. + 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 weights. + Returns an array of the objective targets. """ return [obj.target for obj in self.objectives] From db2ef21029bac25f3875e2680fc84e4d4877fd79 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Sat, 6 Jan 2024 16:29:41 -0500 Subject: [PATCH 15/17] fixed notebooks --- bloptools/bayesian/agent.py | 48 ++++++++++++++++---------- bloptools/bayesian/kernels.py | 40 +++++++++++---------- bloptools/bayesian/models.py | 4 +-- docs/source/tutorials/himmelblau.ipynb | 30 ++++------------ 4 files changed, 59 insertions(+), 63 deletions(-) diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 8c9c08c..5faffa8 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -82,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 @@ -170,7 +170,7 @@ def tell( elif int(self.n_last_trained / self.train_every) != int(len(self.table) / self.train_every): if train_models: - self._train_models(a_priori_hypers=hypers) + self._construct_models(train=train_models, a_priori_hypers=hypers) def _construct_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 @@ -235,11 +235,11 @@ def _construct_models(self, train=True, skew_dims=None, a_priori_hypers=None): # if self.initialized: # self._set_hypers(cached_hypers) # else: - # raise RuntimeError("Could not fit model on initialization!") + # 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): + 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 @@ -261,19 +261,21 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True, upsam 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 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 @@ -332,7 +334,17 @@ def ask(self, acq_func_identifier="qei", n=1, route=True, sequential=True, upsam upsampled_idx = np.linspace(0, len(idx) - 1, upsample * len(idx) - 1) acq_points = sp.interpolate.interp1d(idx, acq_points, axis=0)(upsampled_idx) - return acq_points, acq_func_meta + 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. @@ -366,7 +378,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 @@ -405,8 +417,8 @@ def learn( 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 ---------- @@ -435,9 +447,9 @@ def learn( for i in range(iterations): print(f"running iteration {i + 1} / {iterations}") for single_acq_func in np.atleast_1d(acq_func): - acq_points, acq_func_meta = self.ask(n=n, acq_func_identifier=single_acq_func, upsample=upsample) - new_table = yield from self.acquire(acq_points) - new_table.loc[:, "acq_func"] = acq_func_meta["name"] + 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"] 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} @@ -725,8 +737,8 @@ def all_acq_funcs(self): 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)) 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 fab5f38..a50e84e 100644 --- a/bloptools/bayesian/models.py +++ b/bloptools/bayesian/models.py @@ -15,7 +15,7 @@ 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 ) @@ -31,7 +31,7 @@ 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 ) diff --git a/docs/source/tutorials/himmelblau.ipynb b/docs/source/tutorials/himmelblau.ipynb index 802befa..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(name=\"himmelblau\", target=\"min\")]" + "objectives = [Objective(name=\"himmelblau\", description=\"Himmeblau's function\", target=\"min\")]" ] }, { @@ -170,24 +170,6 @@ "agent.plot_objectives()" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "75a6a54f-0140-4bd0-8367-0dd14cdf8116", - "metadata": {}, - "outputs": [], - "source": [ - "{x % 2 for x in np.arange(10)}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ad31b676-8236-4cf8-b9be-c1214d2e8458", - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "markdown", "id": "dc264346-10fb-4c88-9925-4bfcf0dd3b07", @@ -264,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", ")" ] From 7d82dd14d937f52e67f93bc35d744d8348c2a14a Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Sat, 6 Jan 2024 19:27:30 -0500 Subject: [PATCH 16/17] last PR tweaks --- bloptools/bayesian/agent.py | 288 +++++++++++------- bloptools/bayesian/objectives.py | 7 +- bloptools/bayesian/plotting.py | 6 +- bloptools/bayesian/transforms.py | 2 +- bloptools/tests/conftest.py | 31 ++ bloptools/tests/test_passive_dofs.py | 34 +-- .../{test_agent.py => test_read_write.py} | 0 bloptools/utils/functions.py | 19 ++ ...776a41c3-7a8a-4279-8098-517e1b113fef.ipynb | 6 + 9 files changed, 256 insertions(+), 137 deletions(-) rename bloptools/tests/{test_agent.py => test_read_write.py} (100%) create mode 100644 docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb diff --git a/bloptools/bayesian/agent.py b/bloptools/bayesian/agent.py index 5faffa8..77d1d06 100644 --- a/bloptools/bayesian/agent.py +++ b/bloptools/bayesian/agent.py @@ -126,7 +126,7 @@ def tell( y: Optional[Mapping] = {}, metadata: Optional[Mapping] = {}, append=True, - train_models=True, + train=True, hypers=None, ): """ @@ -145,7 +145,7 @@ def tell( 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. """ if not data: @@ -163,82 +163,153 @@ def tell( self.table = pd.concat([self.table, new_table]) if append else new_table self.table.index = np.arange(len(self.table)) - # TODO: should be a check per model - if len(self.table) > 2: - if any([not hasattr(obj, "model") for obj in self.objectives]): - self._construct_models(train=train_models, a_priori_hypers=hypers) + for obj in self.objectives: + t0 = ttime.monotonic() - elif int(self.n_last_trained / self.train_every) != int(len(self.table) / self.train_every): - if train_models: - self._construct_models(train=train_models, a_priori_hypers=hypers) + cached_hypers = obj.model.state_dict() if hasattr(obj, "model") else None - def _construct_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 + obj.model = self.construct_model(obj) + + 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") + + # TODO: should this be per objective? + self.construct_classifier() + + 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 - inputs = self.table.loc[:, self.dofs.subset(active=True).names].values.astype(float) + def construct_model(self, obj, skew_dims=None): + """ + Construct an untrained model for an objective. + """ - for i, obj in enumerate(self.objectives): - values = self.get_objective_targets(i) - values = np.where(self.all_objectives_valid, values, np.nan) + skew_dims = skew_dims if skew_dims is not None else self.latent_dim_tuples - train_index = ~np.isnan(values) + 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), + ) - if not train_index.sum() >= 2: - raise ValueError("There must be at least two valid data points per objective!") + outcome_transform = botorch.models.transforms.outcome.Standardize(m=1) # , batch_shape=torch.Size((1,))) - train_inputs = torch.tensor(inputs[train_index], dtype=torch.double) - train_values = torch.tensor(values[train_index], dtype=torch.double).unsqueeze(-1) # .unsqueeze(0) + train_inputs = self.train_inputs(active=True) + train_targets = self.train_targets(obj.name) - # for constructing the log normal noise prior - # target_snr = 2e2 - # scale = 2e0 - # loc = np.log(1 / target_snr**2) + scale**2 + safe = ~(torch.isnan(train_inputs).any(axis=1) | torch.isnan(train_targets).any(axis=1)) - 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), - ) + 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, + ) - outcome_transform = botorch.models.transforms.outcome.Standardize(m=1) # , batch_shape=torch.Size((1,))) + model.trained = False - obj.model = models.LatentGP( - train_inputs=train_inputs, - train_targets=train_values, - likelihood=likelihood, - skew_dims=skew_dims, - input_transform=self._subset_input_transform(active=True), - outcome_transform=outcome_transform, - ) + 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 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. @@ -396,22 +467,22 @@ def acquire(self, acquisition_inputs): return products - def load_data(self, data_file, append=True, train_models=True): + 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_models=train_models) + 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_models=True, + 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. @@ -454,7 +525,7 @@ def learn( 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_models=train_models) + 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 @@ -480,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) @@ -517,23 +579,6 @@ def posterior(self, x): def objective_weights_torch(self): return torch.tensor(self.objectives.weights, dtype=torch.double) - def get_objective_targets(self, i): - """Returns the values associated with each objective.""" - - obj = self.objectives[i] - - 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 targets - @property def scalarizing_transform(self): return ScalarizedPosteriorTransform(weights=self.objective_weights_torch, offset=0) @@ -547,27 +592,19 @@ def pseudo_targets(self): """Targets for the posterior transform""" return torch.tensor( [ - self.objectives_targets[..., i].max() + self.train_targets(active=True)[..., i].max() if t == "max" - else self.objectives_targets[..., i].min() + else self.train_targets(active=True)[..., i].min() if t == "min" else t for i, t in enumerate(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(len(self.objectives))], dim=1 - ) - @property def scalarized_objectives(self): """Returns a (n_obs,) array of scalarized objectives""" - return self.targeting_transform.evaluate(self.objectives_targets).sum(axis=-1) - # return (self.objectives_targets * self.objectives.signed_weights).sum(axis=-1) + return self.targeting_transform.evaluate(self.train_targets(active=True)).sum(axis=-1) @property def max_scalarized_objective(self): @@ -632,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) @@ -745,9 +789,47 @@ def all_acq_funcs(self): @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/objectives.py b/bloptools/bayesian/objectives.py index 7bec712..61c8e06 100644 --- a/bloptools/bayesian/objectives.py +++ b/bloptools/bayesian/objectives.py @@ -104,7 +104,12 @@ def __init__(self, objectives: list = []): self.objectives = objectives def __getitem__(self, i): - return self.objectives[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) diff --git a/bloptools/bayesian/plotting.py b/bloptools/bayesian/plotting.py index e51256e..fe5fb9f 100644 --- a/bloptools/bayesian/plotting.py +++ b/bloptools/bayesian/plotting.py @@ -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,7 +87,7 @@ 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): - targets = agent.get_objective_targets(obj_index) + targets = agent.train_targets(obj.name).squeeze(-1).numpy() obj_vmin, obj_vmax = np.nanpercentile(targets, q=[1, 99]) obj_norm = mpl.colors.Normalize(obj_vmin, obj_vmax) diff --git a/bloptools/bayesian/transforms.py b/bloptools/bayesian/transforms.py index 9f64d9c..1b28c6d 100644 --- a/bloptools/bayesian/transforms.py +++ b/bloptools/bayesian/transforms.py @@ -3,7 +3,7 @@ import botorch from botorch.acquisition.objective import PosteriorTransform from botorch.posteriors.gpytorch import GPyTorchPosterior -from botorch.posteriors.posterior_list import PosteriorList # pragma: no cover +from botorch.posteriors.posterior_list import PosteriorList from torch import Tensor diff --git a/bloptools/tests/conftest.py b/bloptools/tests/conftest.py index 034b15c..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 @@ -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 92cfda4..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(name="himmelblau", target="min"), - ] - - 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 100% rename from bloptools/tests/test_agent.py rename to bloptools/tests/test_read_write.py 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-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb b/docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb new file mode 100644 index 0000000..363fcab --- /dev/null +++ b/docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} From 9d0f2d58d07f36de9fed10e356326c8a7e5eb4c1 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Sat, 6 Jan 2024 19:28:13 -0500 Subject: [PATCH 17/17] remove notebook --- ...e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb | 6 ------ 1 file changed, 6 deletions(-) delete mode 100644 docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb diff --git a/docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb b/docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb deleted file mode 100644 index 363fcab..0000000 --- a/docs/source/tutorials/himmelblau-jvsc-fc669c29-d005-463a-8f0e-322008747737776a41c3-7a8a-4279-8098-517e1b113fef.ipynb +++ /dev/null @@ -1,6 +0,0 @@ -{ - "cells": [], - "metadata": {}, - "nbformat": 4, - "nbformat_minor": 5 -}