Multiple non-dominant pareto

[A-ep93]

Hello,
I am trying to redo the tutorial of multi-objective BO posted here. However, when the code reaches to the point where I need to create pareto, I get this error:
BotorchTensorDimensionError: pareto_Y must have the same number of objectives as ref_point. Got 2, expected 1.
So as I understood, the output of pareto_mask = is_non_dominated(train_obj_qparego) should have only 1 true value. However, for me the output of this line can have several true values and I think that is causing this problem. I compared my code with the main code on the website, however, I could not find out what the problem is. I appreciate it if someone can help me to understand what is wrong with my implementation.
This is my code:

import os
import torch
import numpy as np
from botorch.models.gp_regression import SingleTaskGP
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from botorch.models.transforms.outcome import Standardize
from botorch.models import ModelListGP
from gpytorch.mlls.sum_marginal_log_likelihood import SumMarginalLogLikelihood
from botorch import fit_gpytorch_model
from botorch.acquisition.monte_carlo import qExpectedImprovement
from botorch.utils.sampling import sample_simplex
from botorch.acquisition.objective import GenericMCObjective
from botorch.utils.multi_objective.scalarization import get_chebyshev_scalarization
from botorch.sampling.samplers import SobolQMCNormalSampler
from botorch import fit_gpytorch_model
from botorch.acquisition.monte_carlo import qExpectedImprovement, qNoisyExpectedImprovement
from botorch.sampling.samplers import SobolQMCNormalSampler
from botorch.exceptions import BadInitialCandidatesWarning
from botorch.utils.multi_objective.pareto import is_non_dominated
from botorch.utils.multi_objective.hypervolume import Hypervolume
from botorch.optim import optimize_acqf

NUM_RESTARTS = 20
RAW_SAMPLES = 1024

tkwargs = {
    "dtype": torch.double,
    "device": torch.device("cuda" if torch.cuda.is_available() else "cpu"),
}


standard_bounds=torch.tensor([[0,0],[1,1]],dtype=torch.float)


from botorch.test_functions.multi_objective import BraninCurrin
from botorch.utils.sampling import draw_sobol_samples


black_box_function2 = BraninCurrin(negate=True).to(**tkwargs)


def generate_initial_data(n=20):
    train_x = draw_sobol_samples(
        bounds=standard_bounds,n=1, q=n, seed=torch.randint(1000000, (1,)).item()
    ).squeeze(0)
    train_obj = black_box_function2(train_x)
    return train_x, train_obj


def initialize_model(train_x, train_obj):
    models = [
            SingleTaskGP(
                train_x,
                train_obj[..., i : i + 1],
                outcome_transform=Standardize(m=1),
            )
            for i in range(train_obj.shape[-1])
        ]
    model = ModelListGP(*models)
    mll = SumMarginalLogLikelihood(model.likelihood, model)
    return model, mll


def optimize_qparego_and_get_observation(model, train_obj, sampler):
    weights = sample_simplex(2, **tkwargs).squeeze()
    objective = GenericMCObjective(get_chebyshev_scalarization(weights=weights, Y=train_obj))
    acq_func = qExpectedImprovement(  # pyre-ignore: [28]
            model=model,
            objective=objective,
            best_f=objective(train_obj).max(),
            sampler=sampler,
        )
    candidates, _ = optimize_acqf_list(
        acq_function_list=acq_func,
        bounds=standard_bounds,
        num_restarts=NUM_RESTARTS,
        raw_samples=RAW_SAMPLES,  # used for intialization heuristic
        options={"batch_limit": 5, "maxiter": 200},
    )
    new_x =  unnormalize(candidates.detach(), bounds=standard_bounds)
    new_obj = black_box_function2(new_x)
    return new_x, new_obj



NUM_RESTARTS = 20 
MC_SAMPLES=128
RAW_SAMPLES = 1024 
N_TRIALS = 3 
MC_SAMPLES = 128 


ref=torch.tensor([[-18,-6]],dtype=torch.float)
hv = Hypervolume(ref_point=ref)



for trial in range(1, N_TRIALS + 1):
    torch.manual_seed(trial)
    print(f"\nTrial {trial:>2} of {N_TRIALS} ", end="")
    hvs_qparego= []
    train_x_qparego, train_obj_qparego = generate_initial_data(n=6)
    mll_qparego, model_qparego = initialize_model(train_x_qparego, train_obj_qparego)
    pareto_mask = is_non_dominated(train_obj_qparego)
    pareto_y = train_obj_qparego[pareto_mask]
    volume = hv.compute(pareto_y)
    hvs_qparego.append(volume)
    
    fit_gpytorch_model(mll_qparego)
    qparego_sampler = SobolQMCNormalSampler(num_samples=MC_SAMPLES)
    new_x_qparego, new_obj_qparego = optimize_qparego_and_get_observation(
            model_qparego, train_obj_qparego, qparego_sampler
        )
    train_x_qparego = torch.cat([train_x_qparego, new_x_qparego])
    train_obj_qparego = torch.cat([train_obj_qparego, new_obj_qparego])
    for hvs_list, train_obj in zip(
            (hvs_random, hvs_qparego, hvs_qehvi), 
            (train_obj_random, train_obj_qparego, train_obj_qehvi),
        ):
        pareto_mask = is_non_dominated(train_obj)
        pareto_y = train_obj[pareto_mask]
        
        volume = hv.compute(pareto_y)
        hvs_list.append(volume)
    mll_qparego, model_qparego = initialize_model(train_x_qparego, train_obj_qparego)
    hvs_qparego_all.append(hvs_qparego)
    
        
    
    
    

[Balandat]

I think the issue here is that in your code ref_point is a 1 x 2 tensor, while it should be a 2-element single-dimensional tensor.

[A-ep93]

Thanks a lot @Balandat.