import torch


def log_int_prod_gaussians(N1: dict, N2: dict) -> torch.Tensor:
    """Calculate log integral N(x; m1, cov1) * N(x; m2, cov2) dx

    Args:
        N = {'mean': torch.Tensor(k, d), 'cov': torch.Tensor(k, d, d)}

    Returns: k1-by-k2 matrix
    """
    # get number of components
    k1, d = N1['mean'].shape
    k2, d2 = N2['mean'].shape

    assert d == d2

    # expand size
    m1 = torch.unsqueeze(input=N1['mean'], dim=1).expand(k1, k2, d)
    cov1 = torch.unsqueeze(input=N1['cov'], dim=1).expand(k1, k2, d, d)

    m2 = torch.unsqueeze(input=N2['mean'], dim=0).expand(k1, k2, d)
    cov2 = torch.unsqueeze(input=N2['cov'], dim=0).expand(k1, k2, d, d)

    mvn = torch.distributions.multivariate_normal.MultivariateNormal(
        loc=m1 - m2,
        covariance_matrix=cov1 + cov2
    )

    return mvn.log_prob(value=torch.zeros_like(N2['mean'], device=N2['mean'].device))