Oh, sorry about the error. I’m traveling and didn’t have a chance to actually run the code. Gibbs takes a callable for the lengthscale_func, so you’ll need to do something like lengthscale_func = lambda X: lengthscale. For the variance, you can also take a look at ScaledCov in this PR. The cov_func argument would be your Gibbs covariance. The scaling_func argument works like the Gibbs lengthscale argument.
On point 1, you’re multiplying a stationary covariance function by a scaling function, mathematically this is
where \sigma^2(x) is the variance function (a GP in your case), and \odot is the Hadamard (elementwise) product. \sigma^2(x) \sigma^2(x') is an outer product. You may need to constrain it to be non-negative, by using e.g. \exp(\sigma^2) instead.
On point 2, you need to use gp.Latent whenever you don’t have a GP observed with Gaussian noise in the observations. So that means non-Normal observations, or your case where the lengthscale and variance are GPs.