Optimization with descriptor-carrying inputs and engineered Features

Continuous inputs carrying descriptors can be used in combination with engineered features to incorporate more domain knowledge into an optimization. Example could be a solvent of formulation optimization where per continuous input feature additional information about its properties is available.

For examples in the literature have a look at this paper.

Imports

import pandas as pd
import os
import matplotlib.pyplot as plt
import numpy as np

from bofire.data_models.features.api import ContinuousInput, SumFeature, WeightedSumFeature, ContinuousOutput
from bofire.data_models.domain.api import Inputs, EngineeredFeatures, Outputs, Domain, Constraints
from bofire.benchmarks.api import Himmelblau, FormulationWrapper
from bofire.runners.api import run
from bofire.strategies.api import RandomStrategy, SoboStrategy
from bofire.data_models.surrogates.api import SingleTaskGPSurrogate
import bofire.surrogates.api as surrogates
from bofire.data_models.surrogates.api import BotorchSurrogates
from typing import Optional

SMOKE_TEST = os.environ.get("SMOKE_TEST")

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/torch/jit/_script.py:1491: FutureWarning: `torch.jit.script` is deprecated. Please switch to `torch.compile` or `torch.export`.
  warnings.warn(

Problem Setup

For demo purposes, the formulation wrapper is used around the two dimensional Himmelblau Benchmark to wrap it as a formulation problem with 7 features in total. Per original feature, three new features are introduced and one filler feature. The sum of each feature group is used as (normalized) input for the evaluation of the original Himmelblau function. A seventh feature is introduced as filler to always reach a total of 1 to fulfill the formulation/equality constraint. For more information on this setup have a look at the docstring of FormulationWrapper.

bench = FormulationWrapper(benchmark=Himmelblau(), n_features_per_original_feature=3, n_filler_features=1)

The continuos descriptor features have the following descriptor values for the descriptors x_1 and x_2:

pd.concat(
    [
        feat.descriptors.table([feat.key])
        for feat in bench.domain.inputs.get()
        if feat.descriptors is not None
    ]
)

	x_1	x_2
x_1_0	1.0	0.0
x_1_1	1.0	0.0
x_1_2	1.0	0.0
x_2_0	0.0	1.0
x_2_1	0.0	1.0
x_2_2	0.0	1.0

Benchmarks

Random Strategy

As baseline, we run the random strategy on the problem.

def sample(domain):
    strategy = RandomStrategy.make(domain=domain)
    return strategy.ask(10)


def best(domain: Domain, experiments: pd.DataFrame) -> float:
    return experiments.y.min()

random_results = []
for _ in range(15 if not SMOKE_TEST else 1):
    results = run(
        benchmark=bench,
        strategy_factory=lambda domain: RandomStrategy.make(domain=domain),
        n_iterations=40 if not SMOKE_TEST else 2,
        metric=best,
        initial_sampler=sample,
        n_runs=1,
        n_procs=1,
    )
    random_results.append(results[0][0])

  0%|          | 0/2 [00:00<?, ?it/s]Run 0:   0%|          | 0/2 [00:00<?, ?it/s]Run 0:   0%|          | 0/2 [00:00<?, ?it/s, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:00<00:00,  9.94it/s, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:00<00:00,  9.94it/s, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:00<00:00,  9.94it/s, Current Best:=17.142]Run 0: 100%|██████████| 2/2 [00:00<00:00,  9.90it/s, Current Best:=17.142]Run 0: 100%|██████████| 2/2 [00:00<00:00,  9.88it/s, Current Best:=17.142]

Sobo Strategy

Next we run the SoboStrategy on the problem.

def sample(domain):
    return initial_experiments[domain.inputs.get_keys()].copy()

def best(domain: Domain, experiments: pd.DataFrame) -> float:
    return experiments.y.min()

sobo_results = []
for i in range(15 if not SMOKE_TEST else 1):
    initial_experiments = random_results[i][:10].copy()
    results = run(
        benchmark=bench,
        strategy_factory=lambda domain: SoboStrategy.make(domain=domain),
        n_iterations=40 if not SMOKE_TEST else 2,
        metric=best,
        initial_sampler=sample,
        n_runs=1,
        n_procs=1,
    )
    sobo_results.append(results[0][0])

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/bofire/utils/torch_tools.py:51: UserWarning: The given NumPy array is not writable, and PyTorch does not support non-writable tensors. This means writing to this tensor will result in undefined behavior. You may want to copy the array to protect its data or make it writable before converting it to a tensor. This type of warning will be suppressed for the rest of this program. (Triggered internally at /__w/pytorch/pytorch/torch/csrc/utils/tensor_numpy.cpp:213.)
  return torch.from_numpy(np.ascontiguousarray(df.to_numpy())).to(**tkwargs)
  0%|          | 0/2 [00:00<?, ?it/s]Run 0:   0%|          | 0/2 [00:05<?, ?it/s]Run 0:   0%|          | 0/2 [00:05<?, ?it/s, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:05<00:05,  5.11s/it, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:10<00:05,  5.11s/it, Current Best:=17.142]Run 0:  50%|█████     | 1/2 [00:10<00:05,  5.11s/it, Current Best:=14.895]Run 0: 100%|██████████| 2/2 [00:10<00:00,  5.32s/it, Current Best:=14.895]Run 0: 100%|██████████| 2/2 [00:10<00:00,  5.29s/it, Current Best:=14.895]

Now we run the SoboStrategy again, but with an engineered feature as additional input. The WeightedSumFeature computes the sum over the specified descriptors weighted by the values of the involved original features.

def create_strategy(domain):
    surrogate_data = SingleTaskGPSurrogate(
        inputs=domain.inputs,
        outputs=domain.outputs,
        engineered_features=EngineeredFeatures(
            features=[
                WeightedSumFeature(
                    key="WeightedSum",
                    features=[
                        f.key
                        for f in bench.domain.inputs.get()
                        if f.descriptors is not None
                    ],
                    columns=["x_1", "x_2"],
                    )]),
        #kernel=RBFKernel(features=["x_1", "x_2"], lengthscale_prior=HVARFNER_LENGTHSCALE_PRIOR())
    )
    return SoboStrategy.make(domain=domain, surrogate_specs=BotorchSurrogates(surrogates=[surrogate_data]))

sobo_engineered_results_all = []
for i in range(15 if not SMOKE_TEST else 1):
    initial_experiments = random_results[i][:10].copy()
    results = run(
        benchmark=bench,
        strategy_factory=create_strategy,
        n_iterations=40,
        metric=best,
        initial_sampler=sample,
        n_runs=1,
        n_procs=1,
    )
    sobo_engineered_results_all.append(results[0][0])

  0%|          | 0/40 [00:00<?, ?it/s]Run 0:   0%|          | 0/40 [00:05<?, ?it/s]Run 0:   0%|          | 0/40 [00:05<?, ?it/s, Current Best:=17.142]Run 0:   2%|▎         | 1/40 [00:05<03:20,  5.15s/it, Current Best:=17.142]Run 0:   2%|▎         | 1/40 [00:10<03:20,  5.15s/it, Current Best:=17.142]Run 0:   2%|▎         | 1/40 [00:10<03:20,  5.15s/it, Current Best:=17.142]Run 0:   5%|▌         | 2/40 [00:10<03:15,  5.15s/it, Current Best:=17.142]Run 0:   5%|▌         | 2/40 [00:15<03:15,  5.15s/it, Current Best:=17.142]Run 0:   5%|▌         | 2/40 [00:15<03:15,  5.15s/it, Current Best:=17.142]Run 0:   8%|▊         | 3/40 [00:15<03:11,  5.19s/it, Current Best:=17.142]Run 0:   8%|▊         | 3/40 [00:20<03:11,  5.19s/it, Current Best:=17.142]Run 0:   8%|▊         | 3/40 [00:20<03:11,  5.19s/it, Current Best:=2.990] Run 0:  10%|█         | 4/40 [00:20<03:10,  5.30s/it, Current Best:=2.990]Run 0:  10%|█         | 4/40 [00:26<03:10,  5.30s/it, Current Best:=2.990]Run 0:  10%|█         | 4/40 [00:26<03:10,  5.30s/it, Current Best:=2.990]Run 0:  12%|█▎        | 5/40 [00:26<03:06,  5.32s/it, Current Best:=2.990]Run 0:  12%|█▎        | 5/40 [00:31<03:06,  5.32s/it, Current Best:=2.990]Run 0:  12%|█▎        | 5/40 [00:31<03:06,  5.32s/it, Current Best:=0.543]Run 0:  15%|█▌        | 6/40 [00:31<03:01,  5.33s/it, Current Best:=0.543]Run 0:  15%|█▌        | 6/40 [00:37<03:01,  5.33s/it, Current Best:=0.543]Run 0:  15%|█▌        | 6/40 [00:37<03:01,  5.33s/it, Current Best:=0.543]Run 0:  18%|█▊        | 7/40 [00:37<02:57,  5.37s/it, Current Best:=0.543]Run 0:  18%|█▊        | 7/40 [00:42<02:57,  5.37s/it, Current Best:=0.543]Run 0:  18%|█▊        | 7/40 [00:42<02:57,  5.37s/it, Current Best:=0.543]Run 0:  20%|██        | 8/40 [00:42<02:50,  5.32s/it, Current Best:=0.543]Run 0:  20%|██        | 8/40 [00:47<02:50,  5.32s/it, Current Best:=0.543]Run 0:  20%|██        | 8/40 [00:47<02:50,  5.32s/it, Current Best:=0.543]Run 0:  22%|██▎       | 9/40 [00:47<02:44,  5.31s/it, Current Best:=0.543]Run 0:  22%|██▎       | 9/40 [00:53<02:44,  5.31s/it, Current Best:=0.543]Run 0:  22%|██▎       | 9/40 [00:53<02:44,  5.31s/it, Current Best:=0.543]Run 0:  25%|██▌       | 10/40 [00:53<02:42,  5.42s/it, Current Best:=0.543]Run 0:  25%|██▌       | 10/40 [00:58<02:42,  5.42s/it, Current Best:=0.543]Run 0:  25%|██▌       | 10/40 [00:58<02:42,  5.42s/it, Current Best:=0.543]Run 0:  28%|██▊       | 11/40 [00:58<02:38,  5.46s/it, Current Best:=0.543]Run 0:  28%|██▊       | 11/40 [01:03<02:38,  5.46s/it, Current Best:=0.543]Run 0:  28%|██▊       | 11/40 [01:03<02:38,  5.46s/it, Current Best:=0.543]Run 0:  30%|███       | 12/40 [01:03<02:26,  5.24s/it, Current Best:=0.543]Run 0:  30%|███       | 12/40 [01:10<02:26,  5.24s/it, Current Best:=0.543]Run 0:  30%|███       | 12/40 [01:10<02:26,  5.24s/it, Current Best:=0.543]Run 0:  32%|███▎      | 13/40 [01:10<02:37,  5.82s/it, Current Best:=0.543]Run 0:  32%|███▎      | 13/40 [01:17<02:37,  5.82s/it, Current Best:=0.543]Run 0:  32%|███▎      | 13/40 [01:17<02:37,  5.82s/it, Current Best:=0.543]Run 0:  35%|███▌      | 14/40 [01:17<02:34,  5.95s/it, Current Best:=0.543]Run 0:  35%|███▌      | 14/40 [01:22<02:34,  5.95s/it, Current Best:=0.543]Run 0:  35%|███▌      | 14/40 [01:22<02:34,  5.95s/it, Current Best:=0.543]Run 0:  38%|███▊      | 15/40 [01:22<02:25,  5.83s/it, Current Best:=0.543]Run 0:  38%|███▊      | 15/40 [01:28<02:25,  5.83s/it, Current Best:=0.543]Run 0:  38%|███▊      | 15/40 [01:28<02:25,  5.83s/it, Current Best:=0.543]Run 0:  40%|████      | 16/40 [01:28<02:18,  5.77s/it, Current Best:=0.543]Run 0:  40%|████      | 16/40 [01:33<02:18,  5.77s/it, Current Best:=0.543]Run 0:  40%|████      | 16/40 [01:33<02:18,  5.77s/it, Current Best:=0.543]Run 0:  42%|████▎     | 17/40 [01:33<02:11,  5.71s/it, Current Best:=0.543]Run 0:  42%|████▎     | 17/40 [01:40<02:11,  5.71s/it, Current Best:=0.543]Run 0:  42%|████▎     | 17/40 [01:40<02:11,  5.71s/it, Current Best:=0.543]Run 0:  45%|████▌     | 18/40 [01:40<02:10,  5.94s/it, Current Best:=0.543]Run 0:  45%|████▌     | 18/40 [01:46<02:10,  5.94s/it, Current Best:=0.543]Run 0:  45%|████▌     | 18/40 [01:46<02:10,  5.94s/it, Current Best:=0.543]Run 0:  48%|████▊     | 19/40 [01:46<02:03,  5.89s/it, Current Best:=0.543]Run 0:  48%|████▊     | 19/40 [01:52<02:03,  5.89s/it, Current Best:=0.543]Run 0:  48%|████▊     | 19/40 [01:52<02:03,  5.89s/it, Current Best:=0.543]Run 0:  50%|█████     | 20/40 [01:52<01:59,  5.95s/it, Current Best:=0.543]Run 0:  50%|█████     | 20/40 [01:58<01:59,  5.95s/it, Current Best:=0.543]Run 0:  50%|█████     | 20/40 [01:58<01:59,  5.95s/it, Current Best:=0.543]Run 0:  52%|█████▎    | 21/40 [01:58<01:53,  5.99s/it, Current Best:=0.543]Run 0:  52%|█████▎    | 21/40 [02:04<01:53,  5.99s/it, Current Best:=0.543]Run 0:  52%|█████▎    | 21/40 [02:04<01:53,  5.99s/it, Current Best:=0.543]Run 0:  55%|█████▌    | 22/40 [02:04<01:50,  6.12s/it, Current Best:=0.543]Run 0:  55%|█████▌    | 22/40 [02:10<01:50,  6.12s/it, Current Best:=0.543]Run 0:  55%|█████▌    | 22/40 [02:10<01:50,  6.12s/it, Current Best:=0.543]Run 0:  57%|█████▊    | 23/40 [02:10<01:44,  6.12s/it, Current Best:=0.543]Run 0:  57%|█████▊    | 23/40 [02:16<01:44,  6.12s/it, Current Best:=0.543]Run 0:  57%|█████▊    | 23/40 [02:16<01:44,  6.12s/it, Current Best:=0.543]Run 0:  60%|██████    | 24/40 [02:16<01:35,  5.94s/it, Current Best:=0.543]Run 0:  60%|██████    | 24/40 [02:22<01:35,  5.94s/it, Current Best:=0.543]Run 0:  60%|██████    | 24/40 [02:22<01:35,  5.94s/it, Current Best:=0.543]Run 0:  62%|██████▎   | 25/40 [02:22<01:31,  6.12s/it, Current Best:=0.543]Run 0:  62%|██████▎   | 25/40 [02:32<01:31,  6.12s/it, Current Best:=0.543]Run 0:  62%|██████▎   | 25/40 [02:32<01:31,  6.12s/it, Current Best:=0.543]Run 0:  65%|██████▌   | 26/40 [02:32<01:41,  7.28s/it, Current Best:=0.543]Run 0:  65%|██████▌   | 26/40 [02:39<01:41,  7.28s/it, Current Best:=0.543]Run 0:  65%|██████▌   | 26/40 [02:39<01:41,  7.28s/it, Current Best:=0.543]Run 0:  68%|██████▊   | 27/40 [02:39<01:31,  7.04s/it, Current Best:=0.543]Run 0:  68%|██████▊   | 27/40 [02:45<01:31,  7.04s/it, Current Best:=0.543]Run 0:  68%|██████▊   | 27/40 [02:45<01:31,  7.04s/it, Current Best:=0.543]Run 0:  70%|███████   | 28/40 [02:45<01:21,  6.79s/it, Current Best:=0.543]Run 0:  70%|███████   | 28/40 [02:52<01:21,  6.79s/it, Current Best:=0.543]Run 0:  70%|███████   | 28/40 [02:52<01:21,  6.79s/it, Current Best:=0.543]Run 0:  72%|███████▎  | 29/40 [02:52<01:14,  6.76s/it, Current Best:=0.543]Run 0:  72%|███████▎  | 29/40 [02:58<01:14,  6.76s/it, Current Best:=0.543]Run 0:  72%|███████▎  | 29/40 [02:58<01:14,  6.76s/it, Current Best:=0.543]Run 0:  75%|███████▌  | 30/40 [02:58<01:06,  6.68s/it, Current Best:=0.543]Run 0:  75%|███████▌  | 30/40 [03:04<01:06,  6.68s/it, Current Best:=0.543]Run 0:  75%|███████▌  | 30/40 [03:04<01:06,  6.68s/it, Current Best:=0.309]Run 0:  78%|███████▊  | 31/40 [03:04<00:58,  6.46s/it, Current Best:=0.309]Run 0:  78%|███████▊  | 31/40 [03:10<00:58,  6.46s/it, Current Best:=0.309]Run 0:  78%|███████▊  | 31/40 [03:10<00:58,  6.46s/it, Current Best:=0.112]Run 0:  80%|████████  | 32/40 [03:10<00:50,  6.30s/it, Current Best:=0.112]Run 0:  80%|████████  | 32/40 [03:18<00:50,  6.30s/it, Current Best:=0.112]Run 0:  80%|████████  | 32/40 [03:18<00:50,  6.30s/it, Current Best:=0.112]Run 0:  82%|████████▎ | 33/40 [03:18<00:47,  6.78s/it, Current Best:=0.112]Run 0:  82%|████████▎ | 33/40 [03:24<00:47,  6.78s/it, Current Best:=0.112]Run 0:  82%|████████▎ | 33/40 [03:24<00:47,  6.78s/it, Current Best:=0.112]Run 0:  85%|████████▌ | 34/40 [03:24<00:39,  6.56s/it, Current Best:=0.112]Run 0:  85%|████████▌ | 34/40 [03:30<00:39,  6.56s/it, Current Best:=0.112]Run 0:  85%|████████▌ | 34/40 [03:30<00:39,  6.56s/it, Current Best:=0.112]Run 0:  88%|████████▊ | 35/40 [03:30<00:32,  6.42s/it, Current Best:=0.112]Run 0:  88%|████████▊ | 35/40 [03:36<00:32,  6.42s/it, Current Best:=0.112]Run 0:  88%|████████▊ | 35/40 [03:36<00:32,  6.42s/it, Current Best:=0.112]Run 0:  90%|█████████ | 36/40 [03:36<00:25,  6.26s/it, Current Best:=0.112]Run 0:  90%|█████████ | 36/40 [03:43<00:25,  6.26s/it, Current Best:=0.112]Run 0:  90%|█████████ | 36/40 [03:43<00:25,  6.26s/it, Current Best:=0.112]Run 0:  92%|█████████▎| 37/40 [03:43<00:19,  6.39s/it, Current Best:=0.112]Run 0:  92%|█████████▎| 37/40 [03:49<00:19,  6.39s/it, Current Best:=0.112]Run 0:  92%|█████████▎| 37/40 [03:49<00:19,  6.39s/it, Current Best:=0.112]Run 0:  95%|█████████▌| 38/40 [03:49<00:12,  6.24s/it, Current Best:=0.112]Run 0:  95%|█████████▌| 38/40 [03:55<00:12,  6.24s/it, Current Best:=0.112]Run 0:  95%|█████████▌| 38/40 [03:55<00:12,  6.24s/it, Current Best:=0.112]Run 0:  98%|█████████▊| 39/40 [03:55<00:06,  6.22s/it, Current Best:=0.112]Run 0:  98%|█████████▊| 39/40 [04:00<00:06,  6.22s/it, Current Best:=0.112]Run 0:  98%|█████████▊| 39/40 [04:00<00:06,  6.22s/it, Current Best:=0.112]Run 0: 100%|██████████| 40/40 [04:00<00:00,  6.03s/it, Current Best:=0.112]Run 0: 100%|██████████| 40/40 [04:00<00:00,  6.02s/it, Current Best:=0.112]

Plot the performance

Next we plot the performance. Using the engineered features shows a clear benefit. Using a SAAS model on it would lead to even better results, as only two features are needed and SAAS models are very good in figuring out the features of interest.

def get_padded_trajectories(dfs: list[pd.DataFrame]) -> np.ndarray:
    trajectories = []
    for df in dfs:
        y_values = [df.y[:10].max()] + df.y[10:].tolist()
        cummin = np.minimum.accumulate(np.log10(y_values))
        trajectories.append(np.array(cummin))
    max_len = max(len(t) for t in trajectories)
    return [np.pad(traj, (0, max_len-len(traj)), mode='edge') for traj in trajectories]

fig, ax = plt.subplots()

random_trajectories = get_padded_trajectories(random_results)
std_random = np.std(random_trajectories, axis=0)
mean_random = np.mean(random_trajectories, axis=0)
sem_random = std_random / np.sqrt(len(random_trajectories))
ax.plot(mean_random, label='Random', color='blue')
ax.fill_between(range(len(mean_random)), mean_random - sem_random, mean_random + sem_random, color='blue', alpha=0.2)


sobo_trajectories = get_padded_trajectories(sobo_results)
std_sobo = np.std(sobo_trajectories, axis=0)
mean_sobo = np.mean(sobo_trajectories, axis=0)
sem_sobo = std_sobo / np.sqrt(len(sobo_trajectories))
ax.plot(mean_sobo, label='SOBO', color='green')
ax.fill_between(range(len(mean_sobo)), mean_sobo - sem_sobo, mean_sobo + sem_sobo, color='green', alpha=0.2)

sobo_engineered_trajectories = get_padded_trajectories(sobo_engineered_results_all)
std_sobo_engineered = np.std(sobo_engineered_trajectories, axis=0)
mean_sobo_engineered = np.mean(sobo_engineered_trajectories, axis=0)
sem_sobo_engineered = std_sobo_engineered / np.sqrt(len(sobo_engineered_trajectories))
ax.plot(mean_sobo_engineered, label='SOBO Engineered', color='red')
ax.fill_between(range(len(mean_sobo_engineered)), mean_sobo_engineered - sem_sobo_engineered, mean_sobo_engineered + sem_sobo_engineered, color='red', alpha=0.2)

ax.set_title('Optimization Performance')
ax.legend()
plt.show()