random_state not working as expected when using constraints

[IrHaque]

I am using the following versions of the dependencies:
scipy==1.6.2,
scikit-learn==1.4.1.post1,
numpy==1.22.4,
colorama==0.4.6,

on the latest version of BayesianOptimization (1.4.3). The script is run on a cluster set up with HTCondor (I can gather more information regarding the cluster if required).

CONTEXT: I am trying to optimize a function with 5 input parameters $f(\theta_{hS}, \theta_{hX}, \theta_{SX}, v_{S}, v_{X})$. This function has internal constraints, meaning that if a set of input parameters fail the constraints, the function $f$ returns nothing.

I have set up BayesianOptimization using the nonlinear constraint function from scipy scipy.optimize.NonlinearConstraint as described in the section describing advanced constraints in constraints.ipynb.

THE BUG: my issue is that optimizer.max returns slightly different values for each run whilst keeping everything the same between the runs -- I have changed nothing between the runs, I am running the exact same script but getting slightly different results... Below are the outputs I get between the runs.

1.9774822805041e-0 and 1.9741709457416507e-0

As can be seen the latter decimals differs quite a bit, but the overall order of magnitude remains the same.

I have set random_state=1, so that the seed is the same, hence there should be no randomness affecting the optimization. Below I have included my script

#!/usr/bin/python3 python3
import configurer as config
import functions

import os
import subprocess

import numpy as np
import pandas

from scipy.optimize import NonlinearConstraint

if __name__ == '__main__':

    def TRSM(mH1, mH2, mH3, thetahS, thetahX, thetaSX, vs, vx, SM1, SM2, mode, pathTRSM, pathOutput, **kwargs):

        ####################### kwargs #######################

        if 'BFB' in kwargs:
            BFB = kwargs['BFB']

        else: BFB = 0

        if 'Uni' in kwargs:
            Uni = kwargs['Uni']

        else: Uni = 0

        if 'STU' in kwargs:
            STU = kwargs['STU']

        else: STU = 0

        if 'Higgs' in kwargs:
            Higgs = kwargs['Higgs']

        else: Higgs = 0

        ######################################################
        
        # input parameters
        dictModelParams = {'mH1_lb': mH1, 'mH1_ub': mH1,
                           'mH2_lb': mH2, 'mH2_ub': mH2,
                           'mH3_lb': mH3, 'mH3_ub': mH3,
                           'thetahS_lb': thetahS, 'thetahS_ub': thetahS,
                           'thetahX_lb': thetahX, 'thetahX_ub': thetahX,
                           'thetaSX_lb': thetaSX, 'thetaSX_ub': thetaSX,
                           'vs_lb': vs, 'vs_ub': vs,
                           'vx_lb': vx, 'vx_ub': vx}

        # some temporary files needed for the executable below
        dataId = f'mH1{mH1}_mH2{mH2}_mH3{mH3}'
        pathExecutionConfig = os.path.join(pathOutput, f'x_config_{dataId}.tsv')
        pathExecutionOutput = os.path.join(pathOutput, f'x_output_{dataId}.tsv')
        
        # creates contents of pathExecutionConfig
        config.checkCreatorNew(pathExecutionConfig, dictModelParams)

        # for subprocesses.run below
        runTRSM = [pathTRSM, '--BFB', str(BFB), '--Uni', str(Uni), '--STU', str(STU), '--Higgs', str(Higgs), pathExecutionOutput, 'check', pathExecutionConfig]

        # run the executable
        shell_output = subprocess.run(runTRSM, cwd=pathOutput)

        # calculate physical observables from the output from the executable (pathExecutionOutput)
        H1H2, H1H1, H2H2 = functions.ppXNPSM_massfree(pathExecutionOutput, 'mH1', 'mH2', 'mH3',  SM1, SM2,  normalizationSM=1)
        returnValue = (H1H2[mode])[0]
        
        # Sometimes the output is very small values such that it cannot represent them as floats, hence return as 0.
        if np.isnan(returnValue):
            returnValue = 0

        else: pass

        return returnValue

    def constraintTRSM(mH1, mH2, mH3, thetahS, thetahX, thetaSX, vs, vx, SM1, SM2, mode, pathTRSM, pathOutput):
         
        dictModelParams = {'mH1_lb': mH1, 'mH1_ub': mH1,
                           'mH2_lb': mH2, 'mH2_ub': mH2,
                           'mH3_lb': mH3, 'mH3_ub': mH3,
                           'thetahS_lb': thetahS, 'thetahS_ub': thetahS,
                           'thetahX_lb': thetahX, 'thetahX_ub': thetahX,
                           'thetaSX_lb': thetaSX, 'thetaSX_ub': thetaSX,
                           'vs_lb': vs, 'vs_ub': vs,
                           'vx_lb': vx, 'vx_ub': vx}

        dataId = f'mH1{mH1}_mH2{mH2}_mH3{mH3}'
        pathExecutionConfig = os.path.join(pathOutput, f'x_constraint_config_{dataId}.tsv')
        pathExecutionOutput = os.path.join(pathOutput, f'x_constraint_output_{dataId}.tsv')

        config.checkCreatorNew(pathExecutionConfig, dictModelParams)
        runTRSM = [pathTRSM, '--BFB', '1', '--Uni', '1', '--STU', '1', '--Higgs', '1', pathExecutionOutput, 'check', pathExecutionConfig]
        
        shell_output = subprocess.run(runTRSM, cwd=pathOutput)

        try:
            # if returnvalue == 1, then model parameters pass the constraints
            H1H2, H1H1, H2H2 = functions.ppXNPSM_massfree(pathExecutionOutput, 'mH1', 'mH2', 'mH3',  SM1, SM2,  normalizationSM=1)
            returnValue = 1

        # if outputfile (pathExecutionOutput) is empty, the parameters do not pass the constraints
        except pandas.errors.EmptyDataError:
            # if returnValue == -1, then model parameters fail constraints
            returnValue = -1

        return returnValue
        

    # we do not want to optimize mH1, mH2, mH3 and SM1, SM2, pathTRSM, pathOutput are strings required for TRSM
    def wrapInputBayesianOpt(mH1, mH2, mH3, SM1, SM2, mode, pathTRSM, pathOutput):
        def inputBayesianOpt(thetahS, thetahX, thetaSX, vs, vx):
            return TRSM(mH1, mH2, mH3, thetahS, thetahX, thetaSX, vs, vx, SM1, SM2, mode, pathTRSM, pathOutput)
        return inputBayesianOpt

    def wrapInputBayesianOptConstraint(mH1, mH2, mH3, SM1, SM2, mode, pathTRSM, pathOutput):
        def inputBayesianOptConstraint(thetahS, thetahX, thetaSX, vs, vx):
            return constraintTRSM(mH1, mH2, mH3, thetahS, thetahX, thetaSX, vs, vx, SM1, SM2, mode, pathTRSM, pathOutput)
        return inputBayesianOptConstraint


    from bayes_opt import BayesianOptimization
    from bayes_opt import UtilityFunction
    mH1 = 90
    mH2 = 125.09
    mH3 = 500
    SM1, SM2 = 'bb', 'gamgam'
    mode = 4

    black_box_function = wrapInputBayesianOpt(mH1, mH2, mH3, SM1, SM2, mode,
                                              '/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/TRSMBroken', 
                                              '/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt')
    
    black_box_constraint = wrapInputBayesianOptConstraint(mH1, mH2, mH3, SM1, SM2, mode,
                                              '/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/TRSMBroken', 
                                              '/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt')    
    
    # Bounded region of parameter space
    pbounds = {'thetahS': (-np.pi/2, np.pi/2), 'thetahX': (-np.pi/2, np.pi/2), 'thetaSX': (-np.pi/2, np.pi/2), 
               'vs': (1, 1000), 'vx': (1, 1000)}

    constraint = NonlinearConstraint(black_box_constraint, 0, np.inf)

    optimizer = BayesianOptimization(
        f=black_box_function,
        constraint=constraint,
        pbounds=pbounds,
        random_state=1,
        verbose=0,
        allow_duplicate_points=True,
    )

    init_points = 250 
    n_iter = 5
    kind = "ucb"
    kappa='default'
    xi='default'

    print(f'running n_iter = {n_iter}, init_points = {init_points}, kind = {kind}, kappa = {kappa}, xi = {xi}')
    acquisition_function = UtilityFunction(kind=kind)    
    optimizer.maximize(
        init_points=init_points,
        n_iter=n_iter,
        acquisition_function=acquisition_function
    )

    BayOutput = (optimizer.max)

    # IGNORE BELOW, WRITES OUTPUT FROM optimizer.max INTO A .json FILE AND ITERATIONS IN A .tsv FILE

    BayOutput['params']['mH1'] = mH1
    BayOutput['params']['mH2'] = mH2
    BayOutput['params']['mH3'] = mH3
    BayOutput['params']['mode'] = mode

    print(BayOutput)

    BayOutputToJson = {'BayOutput': BayOutput, 'init_points': init_points, 'n_iter': n_iter, 'kind': kind, 'kappa': kappa, 'xi':xi}

    import json
    # save settings used in BayesianOptimization and optimizer.maximize and other useful info
    with open("/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt/x_BayOptOutput.json", "a") as outfile: 
        json.dump(BayOutputToJson, outfile, indent=4)
    
    # clear old contents of iteration file
    with open(f'/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt/x_{mH1}_{mH2}_{mH3}.tsv', 'w') as logs:
        pass
    
    # save iterations
    with open(f'/aaa/bbb/ccc/ddd/eee/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt/x_{mH1}_{mH2}_{mH3}.tsv', 'a') as logs:
        logs.write('iteration\ttarget\tconstraint\tmH1\tmH2\tmH3\tthetahS\tthetahX\tthetaSX\tvs\tvx\n')
        for i, res in enumerate(optimizer.res):
            #logs.write(f"Iteration {i}: \n\t{res}\n")
            logs.write(f'{i}\t{res["target"]}\t{res["constraint"]}\t{BayOutput["params"]["mH1"]}\t{BayOutput["params"]["mH2"]}\t{BayOutput["params"]["mH3"]}\t{res["params"]["thetahS"]}\t{res["params"]["thetahX"]}\t{res["params"]["thetaSX"]}\t{res["params"]["vs"]}\t{res["params"]["vx"]}\n')

I can provide more information about the executable and the modules I am using if needed.

[till-m]

Hey @IrHaque,

I tried to make a more minimal example, but I cannot reproduce your error -- see the snippet below.

Are you sure the function you are optimizing is deterministic?

Assuming that is the case, there could be an issue with our code that my example is too simple to produce. The two major sources of randomness are the GP fitting step and the random sampling step of the acquisition function optimization. You could check that the results produced there are identical.

If you can reduce your example significantly while still retaining the problem, I can also have another look.

import numpy as np
from bayes_opt import BayesianOptimization
import matplotlib.pyplot as plt
from scipy.optimize import NonlinearConstraint

maxs_ = []
for _ in range(10):
    def target_function(x, y):
        # Gardner is looking for the minimum, but this packages looks for maxima, thus the sign switch
        return np.cos(2*x)*np.cos(y) + np.sin(x)

    def constraint_function(x, y):
        return np.cos(x) * np.cos(y) - np.sin(x) * np.sin(y)

    constraint_limit = 0.5

    constraint = NonlinearConstraint(constraint_function, -np.inf, constraint_limit)

    # Bounded region of parameter space
    pbounds = {'x': (0, 6), 'y': (0, 6)}

    optimizer = BayesianOptimization(
        f=target_function,
        constraint=constraint,
        pbounds=pbounds,
        verbose=0, # verbose = 1 prints only when a maximum is observed, verbose = 0 is silent
        random_state=1,
    )

    optimizer.maximize(
        init_points=2,
        n_iter=5,
    )
    maxs_.append(optimizer.max['target'])
print(maxs_) 

# Out: [1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493, 1.7926580410885493

[IrHaque]

Thank you for taking your time checking out my issue @till-m.

I have run your code snippet and I do not see the same issue there. I am fairly sure my function is deterministic, though, I will have to double check that again.

I have been since trying locate the problem and reduce it to a more simple reproducible code snippet as you suggested, however, when I ran my code again with a smaller number of iterations and initial points (eg. init_points = 2, n_iter = 2), I now get the following error message:

Traceback (most recent call last):
  File "/afs/cern.ch/user/i/ihaque/.local/lib/python3.9/site-packages/bayes_opt/bayesian_optimization.py", line 369, in maximize
    x_probe = next(self._queue)
  File "/afs/cern.ch/user/i/ihaque/.local/lib/python3.9/site-packages/bayes_opt/bayesian_optimization.py", line 40, in __next__
    raise StopIteration("Queue is empty, no more objects to retrieve.")
StopIteration: Queue is empty, no more objects to retrieve.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/afs/cern.ch/user/i/ihaque/scannerS/ScannerS-master/build/sh-bbyy-pheno/testing/AtlasLimitsMax_BayesianOpt/tb.py", line 161, in <module>
    optimizer.maximize(
  File "/afs/cern.ch/user/i/ihaque/.local/lib/python3.9/site-packages/bayes_opt/bayesian_optimization.py", line 372, in maximize
    x_probe = self.suggest(util)
  File "/afs/cern.ch/user/i/ihaque/.local/lib/python3.9/site-packages/bayes_opt/bayesian_optimization.py", line 276, in suggest
    y_max_params=self._space.params_to_array(self._space.max()['params']))
TypeError: 'NoneType' object is not subscriptable

I will come back to you if I find what is causing the above.

[till-m]

I now get the following error message:

Quick guess, but it could be that none of the randomly sampled init_points are allowed since they don't fulfill the constraints (thus there is no maximum). Could you try raising the init_points and see if the problem persists?

[till-m]

Hi @IrHaque,

I assume the problem has been fixed. I will close this issue, feel free to re-open if needed.