Multi-dimensional gaussian process regression

Questions
1

I am using a multi dimensional gaussian process. The code is written in python. Note that ‘X_data’ is a list of 10 independent variables and ‘Y_data’ a dependent variable with a non-linear dependency on all of independent variables.
Y_data = fun_non_lin(X_data)

My questions are: 1. is the implementation correct?
2. Are there better implementations?
3. Where can i find necessary documentation for multidimensional gaussian process?
4. Here it is assumed that mean_func is a constant. In reality it’s a complex non-linear function. how can i introduce it?

def Gauss_pymc3(self, method = None):

    with pm.Model() as GP_Model:


        #Start here for Gaussian Process

        GP_const_func   = pm.gp.mean.Constant(np.mean(self.Y_data))

        Sd_gp_len_scale = list(np.max(self.X_data, axis = 1)-
                               np.min(self.X_data, axis = 1))

        GP_len_scale    = []

        for i in range(len(Sd_gp_len_scale)):

            GP_len_scale.append(pm.HalfNormal('GP_len_scale['+str(i)+']',
                                            sigma = Sd_gp_len_scale[i]))

        GP_cov_func     = pm.gp.cov.ExpQuad(len(Sd_gp_len_scale),ls= GP_len_scale)
        GP              = pm.gp.Marginal(mean_func = GP_const_func, cov_func = GP_cov_func)

        noise_l             = pm.Gamma("noise_l", alpha=2, beta=2)
        GP_cov_func_noise   = pm.gp.cov.Exponential(len(Sd_gp_len_scale), noise_l) \
                                        + pm.gp.cov.WhiteNoise(sigma=0.1)

        Y_data_gp       = GP.marginal_likelihood('Y_data_gp',
                                                 X = np.array(self.X_data).T,
                                                 y = self.Y_data,
                                                 noise=GP_cov_func_noise)

        #GP_map          = pm.find_MAP()
        GP_trace        = pm.sample()