Bayesian model calibration with Gaussian Process

Questions
15

Here is a little toy simulation of what I am trying to code:

import numpy as np
import pymc3 as pm

np.random.seed(123)
def black_box(input):
    '''Represents the expensive model f which I am trying to calibrate'''
    l_scale = np.array([1, 1])
    beta = 1
    
    cov_fcn = pm.gp.cov.ExpQuad(2, l_scale)
    cov_mat = cov_fcn(input).eval()
    mean = np.repeat(beta, len(input)) #For simplicity, the mean is constant
    
    return np.random.multivariate_normal(mean, cov_mat).flatten()

x1 = np.random.randn(10)
x2 = np.random.randn(15)
theta_star = np.random.randn(15)
theta_truth = 2

input_1 = np.array([[x_1, theta_truth] for x_1 in x1])
input_2 = np.array([[x2[i], theta_star[i]] for i in range(len(x2))])

#y1 represents observations for known imputs x1 and unknown calibration parameters
y1 = black_box(input_1)
#y2 represents model runs for known x2 and theta_star
y2 = black_box(input_2)

with pm.Model() as toy_model:
    
    #Priors
    theta = pm.Normal("theta", mu = 0, sd = 1)
    beta = pm.Normal("beta", mu = 0, sd = 1)
    l_scale = pm.Gamma("l_scale", alpha= 1, beta = 5, shape = 2)
    σ = pm.HalfCauchy("σ", beta=5)
    
    #Mean and cov of the GP
    mean_fcn = pm.gp.mean.Constant(beta)
    cov_fcn = pm.gp.cov.ExpQuad(2, l_scale)
    
    #This chunk is where it fails###################
    input_1 = np.array([[x_i, theta] for x_i in x1])
    ################################################
    
    gp1 = pm.gp.Marginal(mean_func = mean_fcn, cov_func = cov_fcn)
    gp2 = pm.gp.Marginal(mean_func = mean_fcn, cov_func = cov_fcn)
    
    y_1 = gp2.marginal_likelihood("y1", X=input_1, y=y1, noise=σ)
    y_2 = gp2.marginal_likelihood("y2", X=input_2, y=y2, noise=σ)

This gives me the following error:
AsTensorError: ('Cannot convert [[-1.0856306033005612 theta]\n [0.9973454465835858 theta]\n [0.28297849805199204 theta]\n [-1.506294713918092 theta]\n [-0.5786002519685364 theta]\n [1.651436537097151 theta]\n [-2.426679243393074 theta]\n [-0.42891262885617726 theta]\n [1.265936258705534 theta]\n [-0.8667404022651017 theta]] to TensorType', <class 'numpy.ndarray'>)

Any Idea how to make PyMC3 accept something like:

theta = pm.Normal("theta", mu = 0, sd = 1)
input_1 = np.array([[x_i, theta] for x_i in x1])