How to use the emukit.model_wrappers.gpy_quadrature_wrappers.RBFGPy function in emukit

def qrbf_uniform_infinite():
    x1 = np.array([[-1, 1], [0, 0], [-2, 0.1]])
    x2 = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 3]])
    M1 = x1.shape[0]
    M2 = x2.shape[0]
    D = x1.shape[1]

    # measure
    measure_bounds = [(0, 2), (-4, 3)]
    measure = UniformMeasure(bounds=measure_bounds)

    gpy_kernel = GPy.kern.RBF(input_dim=D)
    emukit_rbf = RBFGPy(gpy_kernel)
    emukit_qrbf = QuadratureRBFUniformMeasure(rbf_kernel=emukit_rbf, integral_bounds=None, measure=measure)
    return emukit_qrbf, x1, x2, M1, M2, D

def qrbf_iso_gauss_infinite():
    x1 = np.array([[-1, 1], [0, 0], [-2, 0.1]])
    x2 = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 3]])
    M1 = x1.shape[0]
    M2 = x2.shape[0]
    D = x1.shape[1]

    gpy_kernel = GPy.kern.RBF(input_dim=D)
    emukit_rbf = RBFGPy(gpy_kernel)
    measure = IsotropicGaussianMeasure(mean=np.arange(D), variance=2.)
    emukit_qrbf = QuadratureRBFIsoGaussMeasure(rbf_kernel=emukit_rbf, measure=measure)
    return emukit_qrbf, x1, x2, M1, M2, D

def qrbf_lebesgue_finite():
    x1 = np.array([[-1, 1], [0, 0], [-2, 0.1]])
    x2 = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 3]])
    M1 = x1.shape[0]
    M2 = x2.shape[0]
    D = x1.shape[1]

    # integral bounds
    bounds = [(-1, 2), (-3, 3)]

    gpy_kernel = GPy.kern.RBF(input_dim=D)
    emukit_rbf = RBFGPy(gpy_kernel)
    emukit_qrbf = QuadratureRBFLebesgueMeasure(emukit_rbf, integral_bounds=bounds)
    return emukit_qrbf, x1, x2, M1, M2, D

def qrbf_uniform_finite():
    x1 = np.array([[-1, 1], [0, 0], [-2, 0.1]])
    x2 = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 3]])
    M1 = x1.shape[0]
    M2 = x2.shape[0]
    D = x1.shape[1]

    # measure
    measure_bounds = [(1, 2), (-4, 2)]
    measure = UniformMeasure(bounds=measure_bounds)

    # integral bounds
    bounds = [(-1, 2), (-3, 3)]

    gpy_kernel = GPy.kern.RBF(input_dim=D)
    emukit_rbf = RBFGPy(gpy_kernel)
    emukit_qrbf = QuadratureRBFUniformMeasure(rbf_kernel=emukit_rbf, integral_bounds=bounds, measure=measure)
    return emukit_qrbf, x1, x2, M1, M2, D

def model_with_density():
    rng = np.random.RandomState(42)
    x_init = rng.rand(5, 2)
    y_init = rng.rand(5, 1)

    gpy_kernel = GPy.kern.RBF(input_dim=x_init.shape[1])
    gpy_model = GPy.models.GPRegression(X=x_init, Y=y_init, kernel=gpy_kernel)
    measure = IsotropicGaussianMeasure(mean=np.arange(x_init.shape[1]), variance=2.)
    qrbf = QuadratureRBFIsoGaussMeasure(RBFGPy(gpy_kernel), measure=measure)
    basegp = BaseGaussianProcessGPy(kern=qrbf, gpy_model=gpy_model)
    model = VanillaBayesianQuadrature(base_gp=basegp, X=x_init, Y=y_init)
    return model

if __name__ == "__main__":
    np.random.seed(0)

    # === Choose MEASURE BELOW ======
    #MEASURE_INTBOUNDS = 'Lebesgue-finite'
    #MEASURE_INTBOUNDS = 'GaussIso-infinite'
    #MEASURE_INTBOUNDS = 'Uniform-infinite'
    MEASURE_INTBOUNDS = 'Uniform-finite'
    # === CHOOSE MEASURE ABOVE ======

    x1 = np.array([[-1, 1], [0, 0], [-2, 0.1]])
    x2 = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 3]])
    D = x1.shape[1]

    gpy_kernel = GPy.kern.RBF(input_dim=D)
    emukit_rbf = RBFGPy(gpy_kernel)

    if MEASURE_INTBOUNDS == 'Lebesgue-finite':
        emukit_qrbf = QuadratureRBFLebesgueMeasure(emukit_rbf, integral_bounds=[(-1, 2), (-3, 3)])
    elif MEASURE_INTBOUNDS == 'GaussIso-infinite':
        measure = IsotropicGaussianMeasure(mean=np.arange(D), variance=2.)
        emukit_qrbf = QuadratureRBFIsoGaussMeasure(rbf_kernel=emukit_rbf, measure=measure)
    elif MEASURE_INTBOUNDS == 'Uniform-infinite':
        measure = UniformMeasure(bounds=[(0, 2), (-4, 3)])
        emukit_qrbf = QuadratureRBFUniformMeasure(emukit_rbf, integral_bounds=None, measure=measure)
    elif MEASURE_INTBOUNDS == 'Uniform-finite':
        measure = UniformMeasure(bounds=[(1, 2), (-4, 2)])
        emukit_qrbf = QuadratureRBFUniformMeasure(emukit_rbf, integral_bounds=[(-1, 2), (-3, 3)], measure=measure)
    else:
        raise ValueError('Measure-integral-bounds combination not defined')

    print()

def loop():
    init_size = 5
    x_init = np.random.rand(init_size, 2)
    y_init = np.random.rand(init_size, 1)
    bounds = [(-1, 1), (0, 1)]

    gpy_model = GPy.models.GPRegression(X=x_init, Y=y_init, kernel=GPy.kern.RBF(input_dim=x_init.shape[1],
                                                                                lengthscale=1., variance=1.))
    emukit_qrbf = QuadratureRBFLebesgueMeasure(RBFGPy(gpy_model.kern), integral_bounds=bounds)
    emukit_model = BaseGaussianProcessGPy(kern=emukit_qrbf, gpy_model=gpy_model)
    emukit_method = VanillaBayesianQuadrature(base_gp=emukit_model, X=x_init, Y=y_init)
    emukit_loop = VanillaBayesianQuadratureLoop(model=emukit_method)
    return emukit_loop, init_size, x_init, y_init

def model():
    rng = np.random.RandomState(42)
    x_init = rng.rand(5, 2)
    y_init = rng.rand(5, 1)

    gpy_kernel = GPy.kern.RBF(input_dim=x_init.shape[1])
    gpy_model = GPy.models.GPRegression(X=x_init, Y=y_init, kernel=gpy_kernel)
    qrbf = QuadratureRBFLebesgueMeasure(RBFGPy(gpy_kernel), integral_bounds=x_init.shape[1] * [(-3, 3)])
    basegp = BaseGaussianProcessGPy(kern=qrbf, gpy_model=gpy_model)
    model = VanillaBayesianQuadrature(base_gp=basegp, X=x_init, Y=y_init)
    return model

def convert_gpy_model_to_emukit_model(gpy_model, integral_bounds: List, integral_name: str = '') \
        -> BaseGaussianProcessGPy:
    """
    Wraps a GPy model and returns an emukit quadrature model

    :param gpy_model: A GPy Gaussian process regression model, GPy.models.GPRegression
    :param integral_bounds: List of D tuples, where D is the dimensionality of the integral and the tuples contain the
        lower and upper bounds of the integral i.e., [(lb_1, ub_1), (lb_2, ub_2), ..., (lb_D, ub_D)]
    :param integral_name: the (variable) name(s) of the integral

    :return: emukit model for quadrature witg GPy backend (IBaseGaussianProcessGPy)
    """

    # wrap standard kernel and get quadrature kernel
    if gpy_model.kern.name is 'rbf':
        standard_kernel_emukit = RBFGPy(gpy_model.kern)
        quadrature_kernel_emukit = QuadratureRBF(standard_kernel_emukit, integral_bounds=integral_bounds,
                                                 integral_name=integral_name)
    else:
        raise NotImplementedError("Only the GPy rbf-kernel is supported right now.")

    # wrap the base-gp model
    return BaseGaussianProcessGPy(kern=quadrature_kernel_emukit, gpy_model=gpy_model)