Multiobjective functionality for Brownie Bee

[dk-teknologisk-mon]

This pull request adds supporting functionality for multiobjective optimization in Brownie Bee. Specifically, Brownie Bee needs to be able to handle plotting of the Pareto front on its own and I have therefore built functions that deliver information about the Pareto front location, its uncertainty in both objective values and a function to select a default point on the front to highlight.

A side-effect of me doing this work is that I discovered several places where functionality related to multiobjective optimization was not reproducible. This was true both for multiobjective opt.ask() and for Pareto front calculations. I have therefore made changes to RNG and seeding in a couple of places to fix this as well.

Below follows a simple script demonstrating the intended use-case that the code supports:

import matplotlib.pyplot as plt
import numpy as np

from ProcessOptimizer import Optimizer
from ProcessOptimizer.plots import get_Brownie_Bee_Pareto
from ProcessOptimizer.model_systems import get_model_system
from ProcessOptimizer.utils.utils import get_Pareto_front_compromise

#%% Set up two simple functions to draw data from
gold_model_system = get_model_system('gold_map')
distance_model_system = get_model_system('distance_map', camp_coordinates=(4,10))

opt = Optimizer(gold_model_system.space, n_initial_points=4, n_objectives=2)

gold_model_system.noise_model.set_seed(40)
distance_model_system.noise_model.set_seed(40)

for i in range(16):
    new_dig_site = opt.ask()
    gold_found = gold_model_system.get_score(new_dig_site)
    distance = distance_model_system.get_score(new_dig_site)
    result_list = opt.tell(new_dig_site, [gold_found, distance])

front_x_data, front_y_data, obj1_error, obj2_error = get_Brownie_Bee_Pareto(opt)

plt.figure(figsize=(10, 10))

# Plot observations
plt.scatter(np.array(opt.yi)[:, 0], np.array(opt.yi)[:, 1], color="k")
# Plot the central points in the Pareto front
plt.plot(front_y_data[:, 0], front_y_data[:, 1], color="k")

# Show bands of uncertainty on each objective
plt.fill_between(
    front_y_data[:, 0], 
    front_y_data[:, 1] - obj2_error, 
    front_y_data[:, 1] + obj2_error, 
    alpha=0.2, 
    color="orange",
)
plt.fill_betweenx(
    front_y_data[:, 1],
    front_y_data[:, 0] - obj1_error,
    front_y_data[:, 0] + obj1_error,
    alpha=0.2,
    color="green",
)
plt.xlabel("Objective 1")
plt.ylabel("Objective 2")

# Get best compromise at present and show it
best_idx = get_Pareto_front_compromise(front_y_data)
plt.scatter(front_y_data[best_idx, 0], front_y_data[best_idx, 1], color="r")

The resulting figure looks like this:

[SRFU-NN]

Nice work on the reproducibility.

I have some comments on the functionality. They are mostly to make sure that we are clear about what we want to do, and that we do that.

[dk-teknologisk-mon]

Updated this pull-request with a third function get_Brownie_Bee_1d_plot(result) that returns the raw data that one needs to build 1D dependency plots in a list of lists. This intent of this function is again to move the construction of plots into the front-end of Brownie Bee, rather than sending up graphics from the back-end.

Among other things this allows us more freedom in how we present plotting information and will also make it easier for us to create overlapping 1D plots for multiobjective optimization mode (instead of the double 1D plots in my mock-up above).

The intent is for this function to be used in conjunction with get_Pareto_front_compromise and get_Brownie_Bee_Pareto in a workflow similar to the following:

# A multiobjective optimizer is fitted

front_x_data, front_y_data, obj1_error, obj2_error = get_Brownie_Bee_Pareto(opt)
best_idx = get_Pareto_front_compromise(front_y_data)

# Some more code that generates the Pareto plot

# Get 1D graphs at the highlighted point on the Pareto front
obj1_1D_data = get_Brownie_Bee_1d_plot(result_list[0], x_eval=front_x_data[best_idx])
obj2_1D_data = get_Brownie_Bee_1d_plot(result_list[1], x_eval=front_x_data[best_idx])

# Some more code that generates overlapping 1D dependency graphs at the highlighted point

[dk-teknologisk-mon]

Also, I have no clue why the Python 3.9 test fails in test_doe.py (for push 059b7b8) as none of my changes touch that part of the package?

This is the specific test that fails:

def test_custom_model(sample_space):
    """Test the get_optimal_DOE function with a custom model."""

    custom_model = "x1 + x2 + x1:x2 + pow(x1, 2)"

    design, factor_names = get_optimal_DOE(
        sample_space, 5, res=5, model=custom_model, seed=42
    )

    design = np.asarray(np.asarray(design[:, :2]), dtype=int)

    expected = np.array([[10, 5], [10, -5], [0, 5], [5, -5], [0, -5]])

    assert design.shape == (5, 2)
    assert np.all(design[:, 0] >= 0) and np.all(design[:, 0] <= 10)
    assert np.all(design[:, 1] >= -5) and np.all(design[:, 1] <= 5)
    assert factor_names == ["x1", "x2"]
    np.testing.assert_array_almost_equal(design, expected)

EDIT: And even more confusingly, the same test now passes...

[SRFU-NN]

Also, I have no clue why the Python 3.9 test fails in test_doe.py (for push 059b7b8) as none of my changes touch that part of the package?

This is the specific test that fails:

def test_custom_model(sample_space):
    """Test the get_optimal_DOE function with a custom model."""

    custom_model = "x1 + x2 + x1:x2 + pow(x1, 2)"

    design, factor_names = get_optimal_DOE(
        sample_space, 5, res=5, model=custom_model, seed=42
    )

    design = np.asarray(np.asarray(design[:, :2]), dtype=int)

    expected = np.array([[10, 5], [10, -5], [0, 5], [5, -5], [0, -5]])

    assert design.shape == (5, 2)
    assert np.all(design[:, 0] >= 0) and np.all(design[:, 0] <= 10)
    assert np.all(design[:, 1] >= -5) and np.all(design[:, 1] <= 5)
    assert factor_names == ["x1", "x2"]
    np.testing.assert_array_almost_equal(design, expected)

EDIT: And even more confusingly, the same test now passes...

That test is finicky, we have had problems with it before. @RuneChristensen-NN can perhaps comment? The failed run is here: https://github.com/novonordisk-research/ProcessOptimizer/actions/runs/19503686383/job/55824051220

[RuneChristensen-NN]

Darn. I thought I already fixed that test. Will do the bad thing for now and just comment out the line that causes failure. Will then revisit at a later point.

[dk-teknologisk-mon]

I also fixed the small bug in plot_objective_1d legend generation, since it was simple and we're touching the nearby code anyway.

closes #357