def METHOD_NAME(p_zx: Distr, p_y1z: Distr, q_z1x: Distr, obs_xy: edic, n_mc: int=0, wlogpi: float=1., repar: bool=True) -> tc.Tensor:
    """ For supervised VAE with structure x <- z -> y.
    Observations are supervised (x,y) pairs.
    For unsupervised observations of x data, use `elbo(p_zx, q_z1x, obs_x)` as VAE z -> x. """
    if n_mc == 0:
        q_y1x_logpval = q_z1x.expect(lambda dc: p_y1z.logp(dc,dc), obs_xy, 0, repar) #, reducefn=tc.logsumexp)
        if hasattr(q_z1x, "entropy"): # No difference for Gaussian
            expc_val = q_z1x.expect(lambda dc: p_zx.logp(dc,dc), obs_xy, 0, repar) + q_z1x.entropy(obs_xy)
        else:
            expc_val = q_z1x.expect(lambda dc: p_zx.logp(dc,dc) - q_z1x.logp(dc,dc), obs_xy, 0, repar)
        return wlogpi * q_y1x_logpval + expc_val
    else:
        q_y1x_pval = q_z1x.expect(lambda dc: p_y1z.logp(dc,dc).exp(), obs_xy, n_mc, repar)
        expc_val = q_z1x.expect(lambda dc: p_y1z.logp(dc,dc).exp() * (p_zx.logp(dc,dc) - q_z1x.logp(dc,dc)),
                obs_xy, n_mc, repar)
        return wlogpi * q_y1x_pval.log() + expc_val / q_y1x_pval
        # q_y1x_logpval = q_z1x.expect(lambda dc: p_y1z.logp(dc,dc), obs_xy, n_mc, repar,
        #         reducefn=tc.logsumexp) - math.log(n_mc)
        # expc_logval = q_z1x.expect(lambda dc: p_y1z.logp(dc,dc) + (p_zx.logp(dc,dc) - q_z1x.logp(dc,dc)).log(),
        #         obs_xy, n_mc, repar, reducefn=tc.logsumexp) - math.log(n_mc)
        # return wlogpi * q_y1x_logpval + (expc_logval - q_y1x_logpval).exp()