code: def _prepare_models() -> Tuple[VGP, VGP, VGP]:
      """
      Prepare models to make sure the coregionalized model with diagonal coregion kernel and
      with fixed lengthscales is equivalent with normal GP regression.
      """
      # 1. Two independent VGPs for two sets of data
      k0 = gpflow.kernels.SquaredExponential()
      set_trainable(k0.lengthscales, False)
      k1 = gpflow.kernels.SquaredExponential()
      set_trainable(k1.lengthscales, False)
      vgp0 = VGP(
          (Datum.X[0], Datum.Y[0]),
          kernel=k0,
          mean_function=Constant(),
          likelihood=gpflow.likelihoods.Gaussian(),
          num_latent_gps=1,
      )
      vgp1 = VGP(
          (Datum.X[1], Datum.Y[1]),
          kernel=k1,
          mean_function=Constant(),
          likelihood=gpflow.likelihoods.Gaussian(),
          num_latent_gps=1,
      )
      # 2. Coregionalized VGP
      kc = gpflow.kernels.SquaredExponential(active_dims=[0, 1])
      set_trainable(kc.lengthscales, False)
      set_trainable(
          kc.variance, False
      )  # variance is handled by the Coregion kernel
      coreg = gpflow.kernels.Coregion(
          output_dim=2, rank=1, active_dims=[2]
      )
      coreg.W.assign(
          np.zeros((2, 1))
      )  # zero correlation between outputs
      set_trainable(coreg.W, False)
      lik = gpflow.likelihoods.SwitchedLikelihood(
          [gpflow.likelihoods.Gaussian(), gpflow.likelihoods.Gaussian()]
      )
      mean_c = gpflow.mean_functions.SwitchedMeanFunction(
          [
              gpflow.mean_functions.Constant(),
              gpflow.mean_functions.Constant(),
          ]
      )
      cvgp = VGP(
          (Datum.X_augmented, Datum.Y_augmented),
          kernel=kc * coreg,
          mean_function=mean_c,
          likelihood=lik,
          num_latent_gps=1,
      )
  
      # Train them for a small number of iterations
  
      opt = gpflow.optimizers.Scipy()
      opt.minimize(
          vgp0.training_loss,
          variables=vgp0.trainable_variables,
          options=dict(maxiter=1000),
          method="BFGS",
      )
      opt.minimize(
          vgp1.training_loss,
          variables=vgp1.trainable_variables,
          options=dict(maxiter=1000),
          method="BFGS",
      )
      opt.minimize(
          cvgp.training_loss,
          variables=cvgp.trainable_variables,
          options=dict(maxiter=1000),
          method="BFGS",
      )
  
      return vgp0, vgp1, cvgp
  
file_path: GPflow/tests/gpflow/kernels/test_coregion.py
imports: from typing import Sequence, Tuple
  import numpy
  from check_shapes import ShapeChecker
  from numpy.testing import assert_allclose
  import gpflow
  from gpflow import set_trainable
  from gpflow.base import AnyNDArray
  from gpflow.mean_functions import Constant
  from gpflow.models import VGP
label_desc: Use image node_id 1 for calling a global function with example usage: _prepare_models() and returns: vgp0,  vgp1,  cvgp