def l2(Ks, dim, X_rhos, Y_rhos, required, clamp=True, to_self=False):
n_X = len(X_rhos)<EOL>n_Y = len(Y_rhos)<EOL>linears = required<EOL>assert linears.shape == (<NUM_LIT:1>, Ks.size, n_X, n_Y, <NUM_LIT:2>)<EOL>X_quadratics = np.empty((Ks.size, n_X), dtype=np.float32)<EOL>for i, rho in enumerate(X_rhos):<EOL><INDENT>X_quadratics[:, i] = quadratic(Ks, dim, rho)<EOL><DEDENT>Y_quadratics = np.empty((Ks.size, n_Y), dtype=np.float32)<EOL>for j, rho in enumerate(Y_rhos):<EOL><INDENT>Y_quadratics[:, j] = quadratic(Ks, dim, rho)<EOL><DEDENT>est = -linears.sum(axis=<NUM_LIT:4>)<EOL>est += X_quadratics[None, :, :, None]<EOL>est += Y_quadratics[None, :, None, :]<EOL>np.maximum(est, <NUM_LIT:0>, out=est)<EOL>np.sqrt(est, out=est)<EOL>if to_self:<EOL><INDENT>est[:, :, xrange(n_X), xrange(n_Y)] = <NUM_LIT:0><EOL><DEDENT>return est[:, :, :, :, None]<EOL>
r'''
    Estimates the L2 distance between distributions, via
        \int (p - q)^2 = \int p^2 - \int p q - \int q p + \int q^2.

    \int pq and \int qp are estimated with the linear function (in both
    directions), while \int p^2 and \int q^2 are estimated via the quadratic
    function below.

    Always clamps negative estimates of l2^2 to 0, because otherwise the sqrt
    would break.