[Bug] Penicillin test function gives better values than in the reference paper

[mfeurer]

🐛 Bug

The penicillin test function gives better values than in the reference paper. In the original paper, the authors sample 10**6 configurations and obtain an optimum around -14.65. Doing the same with the BoTorch implementation leads to values smaller than -20.

To reproduce

** Code snippet to reproduce **

import numpy as np
import torch

from botorch.test_functions.multi_objective import Penicillin
from botorch.utils.transforms import unnormalize

p = Penicillin(noise_std=0, negate=True)

# Draw a lot of samples as in the original paper
input_orig = torch.tensor(np.random.random((1000000, 7))).double()
input = unnormalize(input_orig, bounds=torch.tensor(p._bounds).transpose(0, 1))
values = p.evaluate_true(input)[:, 0]
argmin = np.argmin(values.numpy())
print(values[argmin], input[argmin], input_orig[argmin])

# Check whether this value might be due to scalarization
print(p.evaluate_true(input[argmin]))

** Stack trace/error message **

tensor(-23.1735, dtype=torch.float64) tensor([8.6234e+01, 1.1553e+00, 3.0218e+02, 1.7137e+01, 3.2586e-02, 6.2063e+02,
        5.3257e+00], dtype=torch.float64) tensor([0.4372, 0.0616, 0.9176, 0.9519, 0.0461, 0.6032, 0.2171],
       dtype=torch.float64)
tensor([-23.1735,  29.2081, 496.0000], dtype=torch.float64)

Expected Behavior

A value around -14.65

System information

Please complete the following information:

• 0.7.0
• 1.9.0
• 1.13.1+cu117
• Ubuntu 20.04

Additional context

I tried to feed the found value into the original code from https://github.com/HarryQL/TuRBO-Penicillin but that fails with the best value found with the BoTorch Penicillin function:

import numpy as np
import torch

# X: biomass concentration
# g/L
# X_0 = 0.1

# S:substrate concentration
# g/L
# S_0 = 0.1

# CL:dissolved oxygen concentration

# P: penicillin concentration
# g/L
# P_0 = 0


# CO2:carbon dioxide concentration;

# H:hydrogen ion concentration for pH

# T: temperature.

C_L_star = 8.26

Y_xs = 0.45
Y_xo = 0.04
Y_ps = 0.90
Y_po = 0.20

K_1 = 10**(-10)
K_2 = 7 * 10**(-5)
m_X = 0.014
m_o = 0.467

alpha_1 = 0.143
alpha_2 = 4*10**(-7)
alpha_3 = 10**(-4)
mu_X = 0.092
K_X = 0.15
# K_ox = 2*10**(-2)
# K_op = 5*10**(-4)
mu_p = 0.005
K_p = 0.0002
K_I = 0.10
p = 3
K = 0.04
k_g = 7 * 10**(3)
E_g = 5100
k_d = 10**(33)
E_d = 50000

# rou_dot_C_p = 1/1500
# rou_c_dot_C_pc = 1/2000

rou_dot_C_p = 1000
rou_c_dot_C_pc = 1000

r_q1 = 60
r_q2 = 1.6783 * 10**(-4)
a = 1000
b = 0.60

alpha = 70
beta = 0.4
lambd = 2.5 * 10**(-4)
gamma = 10**(-5)

# kelvin
T_v = 273
T_o = 373

# CAL/(MOL K)
R = 1.9872

# total unit time: hrs
t = 2500
V_max = 180
V_limits = [60, 120]
X_limits = [0.05, 18]
CO2 = 0
T_limits = [293, 303]
S_limits = [0.05, 18]
F_limits = [0.01, 0.50]
s_f_limits = [500, 700]
H_limits = [5, 6.5]
limits = [V_limits, X_limits, T_limits, S_limits, F_limits, s_f_limits, H_limits]


def penicilin_exp_BO(X_input):
    #     V_limits, X_limits, T_limits, S_limits, F_limits, s_f_limits, H_limits

    V, X, T, S, F, s_f, H_ = X_input[:, 0], X_input[:, 1], X_input[:, 2], X_input[:, 3], X_input[:, 4], X_input[:, 5], X_input[:, 6]

    P = 0
    CO2 = 0
    t = 2500

    l_P = []
    l_V = []
    l_X = []
    l_T = []
    l_S = []
    l_F = []
    l_s_f = []
    l_H_ = []
    l_CO2 = []
    l_t = []

    l_P.append(P)
    l_V.append(V)
    l_X.append(X)
    l_T.append(T)
    l_S.append(S)
    l_F.append(F)
    l_s_f.append(s_f)
    l_H_.append(H_)
    l_CO2.append(CO2)
    l_t.append(0)

    H = 10 ** (-H_)

    for i in np.arange(t) + 1:

        F_loss = V * lambd * (np.exp(5 * ((T - T_o) / (T_v - T_o))) - 1)
        dV_dt = F - F_loss

        mu = (mu_X / (1 + K_1 / H + H / K_2)) * (S / (K_X * X + S)) * (
                    (k_g * np.exp(-E_g / (R * T))) - (k_d * np.exp(-E_d / (R * T))))
        dX_dt = mu * X - (X / V) * dV_dt

        mu_pp = mu_p * (S / (K_p + S + S ** 2 / K_I))
        dS_dt = - (mu / Y_xs) * X - (mu_pp / Y_ps) * X - m_X * X + F * s_f / V - (S / V) * dV_dt

        dP_dt = (mu_pp * X) - K * P - (P / V) * dV_dt

        dCO2_dt = alpha_1 * dX_dt + alpha_2 * X + alpha_3

        # UPDATE
        P = P + dP_dt
        V = V + dV_dt
        X = X + dX_dt
        S = S + dS_dt
        CO2 = CO2 + dCO2_dt

        l_P.append(P)
        l_V.append(V)
        l_X.append(X)
        l_T.append(T)
        l_S.append(S)
        l_F.append(F)
        l_s_f.append(s_f)
        l_H_.append(H_)
        l_CO2.append(CO2)
        l_t.append(i)

        if V > V_max:
            #             print('Too large V')
            break

        if S < 0:
            #             print('Too small S')
            break

        if dP_dt < 10e-12:
            #             print('Converged P')
            break

    print(P)
    return -P


query = [8.6234e+01, 1.1553e+00, 3.0218e+02, 1.7137e+01, 3.2586e-02, 6.2063e+02, 5.3257e+00]
print(query)
print(
    penicilin_exp_BO(
        torch.tensor(query).reshape((1, -1)).float()
    )
)

which gives

[86.234, 1.1553, 302.18, 17.137, 0.032586, 620.63, 5.3257]
Traceback (most recent call last):
  File "/home/feurerm/projects/2022_pfns_4_bo/local_harryql.py", line 177, in <module>
    penicilin_exp_BO(
  File "/home/feurerm/projects/2022_pfns_4_bo/local_harryql.py", line 130, in penicilin_exp_BO
    (k_g * np.exp(-E_g / (R * T))) - (k_d * np.exp(-E_d / (R * T))))
RuntimeError: Overflow when unpacking long

[Balandat]

@sdaulton any idea what could be going on here? Given the error and looking through the code it may well be that there are some numerical issues going on (in either implementation).