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 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)