Skip to content

Commit

Permalink
implemented qRollout
Browse files Browse the repository at this point in the history
  • Loading branch information
FilippoAiraldi committed Dec 31, 2023
1 parent e87fd3a commit 16b374b
Show file tree
Hide file tree
Showing 4 changed files with 167 additions and 5 deletions.
55 changes: 52 additions & 3 deletions benchmarking/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
import numpy as np
import torch
from botorch.acquisition import ExpectedImprovement
from botorch.acquisition.multi_step_lookahead import warmstart_multistep
from botorch.fit import fit_gpytorch_mll
from botorch.models import SingleTaskGP
from botorch.models.model import Model
Expand All @@ -28,11 +29,15 @@
from status import filter_tasks_by_status
from torch import Tensor

from globopt.myopic_acquisitions import (
from globopt import (
GaussHermiteSampler,
IdwAcquisitionFunction,
make_idw_acq_factory,
qIdwAcquisitionFunction,
qRollout,
)

# from globopt.nonmyopic_acquisitions import qRollout
from globopt.problems import get_available_benchmark_problems, get_benchmark_problem
from globopt.regression import Idw, Rbf

Expand Down Expand Up @@ -169,7 +174,51 @@ def next_obs(model: Union[Rbf, Idw], *_) -> tuple[Tensor, float, Tensor]:
return X_opt, acq_opt.item(), torch.nan

elif method.startswith("rollout"):
raise NotImplementedError
horizon = int(method.split(".")[1])
n_restarts_ = int(n_restarts * (horizon + 1) / 2)
raw_samples_ = int(raw_samples * (horizon + 1) / 2)
maxfun = 15_000
gh_sampler = GaussHermiteSampler(sample_shape=torch.Size([16]))
kwargs_factory = make_idw_acq_factory(c1, c2)

def next_obs(
model: Union[Rbf, Idw], budget: int, prev_full_opt: Tensor
) -> tuple[Tensor, float]:
remaining_horizon = min(horizon, budget)
if remaining_horizon == 1:
acqfun = qIdwAcquisitionFunction(model, c1, c2, sampler=gh_sampler)
X_opt, acq_opt = optimize_acqf(
acqfun, bounds, 1, n_restarts, raw_samples, {"seed": mk_seed()}
)
return X_opt, acq_opt.item(), torch.nan

acqfun = qRollout(
model,
horizon,
qIdwAcquisitionFunction,
kwargs_factory,
valfunc_sampler=gh_sampler,
)
if prev_full_opt is torch.nan:
prev_full_opt = None
else:
prev_full_opt = warmstart_multistep(
acqfun, bounds, n_restarts_, raw_samples_, prev_full_opt
)
full_opt, tree_vals = optimize_acqf(
acqfun,
bounds,
horizon, # == acqfun.get_augmented_q_batch_size(1)
n_restarts_,
raw_samples_,
batch_initial_conditions=prev_full_opt,
return_best_only=False,
return_full_tree=True,
options={"seed": mk_seed(), "maxfun": maxfun},
)
best_tree_idx = tree_vals.argmax()
X_opt = acqfun.extract_candidates(full_opt[best_tree_idx])
return X_opt, tree_vals[best_tree_idx].item(), full_opt

else:
raise NotImplementedError
Expand Down Expand Up @@ -264,7 +313,7 @@ def run_benchmarks(
" `myopic`. Non-myopic algorithms are `rollout` and `multi-tree`, where the "
"horizons to simulate can be specified with a dot folloewd by (one or more) "
"horizons, e.g., `rollout.2.3`. These horizons should be larger than 1. "
"Methods are: `random`, `ei`, `myopic`, `s-myopic`.",
"Methods are: `random`, `ei`, `myopic`, `s-myopic`, `rollout`.",
required=True,
)
group.add_argument(
Expand Down
3 changes: 3 additions & 0 deletions src/globopt/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@

__all__ = [
"qIdwAcquisitionFunction",
"qRollout",
"make_idw_acq_factory",
"GaussHermiteSampler",
"Idw",
"IdwAcquisitionFunction",
Expand All @@ -10,5 +12,6 @@
]

from globopt.myopic_acquisitions import IdwAcquisitionFunction, qIdwAcquisitionFunction
from globopt.nonmyopic_acquisitions import make_idw_acq_factory, qRollout
from globopt.regression import Idw, Rbf
from globopt.sampling import GaussHermiteSampler, PosteriorMeanSampler
8 changes: 6 additions & 2 deletions src/globopt/myopic_acquisitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,9 @@ def __init__(
contribution is null).
"""
super().__init__(model)
Y_min, Y_max = model.train_Y.aminmax(dim=-2, keepdim=True)
# Y_min, Y_max = model.train_Y.aminmax(dim=-2, keepdim=True)
Y_min = model.train_Y.amin(dim=-2, keepdim=True)
Y_max = model.train_Y.amax(dim=-2, keepdim=True)
self.register_buffer("span_Y", (Y_max - Y_min).clamp_min(span_Y_min))
self.register_buffer("c1", torch.scalar_tensor(c1))
self.register_buffer("c2", torch.scalar_tensor(c2))
Expand Down Expand Up @@ -195,7 +197,9 @@ def __init__(
contribution is null).
"""
super().__init__(model, sampler)
Y_min, Y_max = model.train_Y.aminmax(dim=-2, keepdim=True)
# Y_min, Y_max = model.train_Y.aminmax(dim=-2, keepdim=True)
Y_min = model.train_Y.amin(dim=-2, keepdim=True)
Y_max = model.train_Y.amax(dim=-2, keepdim=True)
self.register_buffer("span_Y", (Y_max - Y_min).clamp_min(span_Y_min))
self.register_buffer("c1", torch.scalar_tensor(c1))
self.register_buffer("c2", torch.scalar_tensor(c2))
Expand Down
106 changes: 106 additions & 0 deletions src/globopt/nonmyopic_acquisitions.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,106 @@
from typing import Any, Optional

import torch
from botorch.acquisition import AcquisitionFunction
from botorch.acquisition.multi_step_lookahead import (
TAcqfArgConstructor,
qMultiStepLookahead,
)
from botorch.models.model import Model
from botorch.sampling.base import MCSampler

from globopt.sampling import PosteriorMeanSampler

# per se, this should support `q > 1`; however, it behaves like an analyical
# acquisition function because in practice:
# - we fix the batch sizes to 1, so the `q`'s in the horizon are fix to 1
# - we couple `qIdwAcquisitionFunction` with `PosteriorMeanSampler`, meaning the
# number of fantasies does not exponentially grow


def make_idw_acq_factory(
c1: float, c2: float, span_Y_min: float = 1e-3
) -> TAcqfArgConstructor:
"""Returns a kwargs factory for `IdwAcquisitionFunction` with the given parameters,
useful for `qMultiStepLookahead`."""

def _inner(*_, **__) -> dict[str, Any]:
return {"c1": c1, "c2": c2, "span_Y_min": span_Y_min}

return _inner


class qRollout(qMultiStepLookahead):
"""Rollout nonmyopic acquisition function based on `qMultiStepLookahead`.
This acquisition function rolls out the base acquisition function along the given
horizon, and returns the sum of the values of the base acquisition function at each
stage. Rollout is known to always outperform greedy selection, and is a good tool
for improving the performance of myopic base acquisition functions.
"""

def __init__(
self,
model: Model,
horizon: int,
valfunc_cls: type[AcquisitionFunction], # base policy
valfunc_argfactory: Optional[TAcqfArgConstructor] = None,
batch_sizes: Optional[list[int]] = None, # value of `q` at each stage
fantasies_sampler: Optional[MCSampler] = None,
valfunc_sampler: Optional[MCSampler] = None,
) -> None:
"""Instantiates the rollout acquisition function.
Parameters
----------
model : Model
A fitted model.
horizon : int
Length of the rollout horizon. Must be at least 2.
valfunc_cls : type[AcquisitionFunction]
The type of the base acquisition function class.
valfunc_argfactory: TAcqfArgConstructor, optional
A callable that takes the current model and observatiosn and returns
the kwargs to pass to the base acquisition function constructor.
batch_sizes : list[int], optional
A list `[q_1, ..., q_k]` containing the batch sizes for the `k` look-ahead
steps. By default, all batch sizes are set to 1 along the horizon.
fantasies_sampler : MCSampler, optional
Sampler used to sample the fantasies. By default, `PosterionMeanSampler` is
used, i.e., a sampler that always takes the posterior mean as the single
sample.
valfunc_sampler : MCSampler, optional
A custom sampler to override the sampling of the base acquisition function
values.
"""
# prepare and check inputs
if horizon < 2:
raise ValueError("horizon must be at least 2")
if batch_sizes is None:
batch_sizes = [1] * (horizon - 1)
if fantasies_sampler is None: # by default, sample the posterior mean
fantasies_sampler = PosteriorMeanSampler()
no_valfunc_sampler = valfunc_sampler is None
if no_valfunc_sampler:
inner_mc_sample = None
else:
if len(valfunc_sampler.sample_shape) > 1:
raise ValueError("`valfunc_sampler` must have a single sample shape")
inner_mc_sample = valfunc_sampler.sample_shape[0]

# construct base
super().__init__(
model=model,
batch_sizes=batch_sizes,
samplers=[fantasies_sampler] * (horizon - 1),
valfunc_cls=[valfunc_cls] * horizon,
valfunc_argfacs=[valfunc_argfactory] * horizon,
inner_mc_samples=[inner_mc_sample] * horizon,
)

# override inner samplers post-construction
if not no_valfunc_sampler:
new_samplers = []
for sampler in self.inner_samplers:
new_samplers.append(None if sampler is None else valfunc_sampler)
self.inner_samplers = torch.nn.ModuleList(new_samplers)

0 comments on commit 16b374b

Please sign in to comment.