docs: def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool=False, **Ys):
      """
      Computes N Gaussian expectation integrals of one or more functions
      using Gauss-Hermite quadrature. The Gaussians must be independent.
  
      :param funcs: the integrand(s):
          Callable or Iterable of Callables that operates elementwise, on the following arguments:
          - `Din` positional arguments to match Fmu and Fvar; i.e., 1 if Fmu and Fvar are tensors;
            otherwise len(Fmu) (== len(Fvar)) positional arguments F1, F2, ...
          - the same keyword arguments as given by **Ys
          All arguments will be tensors of shape (N, 1)
  
      :param H: number of Gauss-Hermite quadrature points
      :param Fmu: array/tensor or `Din`-tuple/list thereof
      :param Fvar: array/tensor or `Din`-tuple/list thereof
      :param logspace: if True, funcs are the log-integrands and this calculates
          the log-expectation of exp(funcs)
      :param **Ys: arrays/tensors; deterministic arguments to be passed by name
  
      Fmu, Fvar, Ys should all have same shape, with overall size `N` (i.e., shape (N,) or (N, 1))
      :return: shape is the same as that of the first Fmu
  
  
      Example use-cases:
  
          Fmu, Fvar are mean and variance of the latent GP, can be shape (N, 1) or (N,)
          m1, m2 are 'scalar' functions of a single argument F, broadcasting over arrays
  
          Em1, Em2 = ndiagquad([m1, m2], 50, Fmu, Fvar)
              calculates Em1 = ∫ m1(F) N(F; Fmu, Fvar) dF and Em2 = ∫ m2(F) N(F; Fmu, Fvar) dF
              for each of the elements of Fmu and Fvar. Em1 and Em2 have the same shape as Fmu.
  
          logp is a 'scalar' function of F and Y
          Y are the observations, with shape (N,) or (N, 1) with same length as Fmu and Fvar
  
          Ev = ndiagquad(logp, 50, Fmu, Fvar, Y=Y)
              calculates Ev = ∫ logp(F, Y) N(F; Fmu, Fvar) dF (variational expectations)
              for each of the elements of Y, Fmu and Fvar. Ev has the same shape as Fmu.
  
          Ep = ndiagquad(logp, 50, Fmu, Fvar, logspace=True, Y=Y)
              calculates Ep = log ∫ exp(logp(F, Y)) N(F; Fmu, Fvar) dF (predictive density)
              for each of the elements of Y, Fmu and Fvar. Ep has the same shape as Fmu.
  
          Heteroskedastic likelihoods:
          g1, g2 are now functions of both F and G
          logp is a function of F, G and Y
          Gmu, Gvar are mean and variance of a different GP controlling the variance
  
          Em = ndiagquad(m1, 50, Fmu, Fvar)
              -> Em1 = ∫∫ m1(F, G) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG
  
          Ev = ndiagquad(logp, 50, Fmu, Fvar, Y=Y)
              -> Ev = ∫∫ logp(F, G, Y) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG
                 (variational expectations)
  
          Ep = ndiagquad(logp, 50, Fmu, Fvar, logspace=True, Y=Y)
              -> Ep = log ∫∫ exp(logp(F, G, Y)) N(F; Fmu, Fvar) N(G; Gmu, Gvar) dF dG
                 (predictive density)
      """
      def unify(f_list):
          """
          Stack a list of means/vars into a full block
          """
          return tf.reshape(
                  tf.concat([tf.reshape(f, (-1, 1)) for f in f_list], axis=1),
                  (-1, 1, Din))
  
      if isinstance(Fmu, (tuple, list)):
          Din = len(Fmu)
          shape = tf.shape(Fmu[0])
          Fmu, Fvar = map(unify, [Fmu, Fvar])    # both N x 1 x Din
      else:
          Din = 1
          shape = tf.shape(Fmu)
          Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]]
  
      xn, wn = mvhermgauss(H, Din)
      # xn: H**Din x Din, wn: H**Din
  
      gh_x = xn.reshape(1, -1, Din)             # 1 x H**Din x Din
      Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu   # N x H**Din x Din
      Xs = [Xall[:, :, i] for i in range(Din)]  # N x H**Din  each
  
      gh_w = wn * np.pi ** (-0.5 * Din)  # H**Din x 1
  
      for name, Y in Ys.items():
          Y = tf.reshape(Y, (-1, 1))
          Y = tf.tile(Y, [1, H**Din])  # broadcast Y to match X
          # without the tiling, some calls such as tf.where() (in bernoulli) fail
          Ys[name] = Y  # now N x H**Din
  
      def eval_func(f):
          feval = f(*Xs, **Ys)  # f should be elementwise: return shape N x H**Din
          if logspace:
              log_gh_w = np.log(gh_w.reshape(1, -1))
              result = tf.reduce_logsumexp(feval + log_gh_w, axis=1)
          else:
              result = tf.matmul(feval, gh_w.reshape(-1, 1))
          return tf.reshape(result, shape)
  
      if isinstance(funcs, Iterable):
          return [eval_func(f) for f in funcs]
      else:
          return eval_func(funcs)