def jensen_shannon(Ks, dim, X_rhos, Y_rhos, required,<EOL>clamp=True, to_self=False):
X_ns = np.array([rho.shape[<NUM_LIT:0>] for rho in X_rhos])<EOL>Y_ns = np.array([rho.shape[<NUM_LIT:0>] for rho in Y_rhos])<EOL>n_X = X_ns.size<EOL>n_Y = Y_ns.size<EOL>cores = required<EOL>assert cores.shape == (<NUM_LIT:1>, Ks.size, n_X, n_Y, <NUM_LIT:2>)<EOL>X_bits = np.empty((Ks.size, n_X), dtype=np.float32)<EOL>for i, rho in enumerate(X_rhos):<EOL><INDENT>X_bits[:, i] = dim * np.mean(np.log(rho), axis=<NUM_LIT:0>)<EOL><DEDENT>X_bits += np.log(X_ns - <NUM_LIT:1>)[np.newaxis, :]<EOL>Y_bits = np.empty((Ks.size, n_Y), dtype=np.float32)<EOL>for j, rho in enumerate(Y_rhos):<EOL><INDENT>Y_bits[:, j] = dim * np.mean(np.log(rho), axis=<NUM_LIT:0>)<EOL><DEDENT>Y_bits += np.log(Y_ns - <NUM_LIT:1>)[np.newaxis, :]<EOL>est = cores.sum(axis=<NUM_LIT:4>)<EOL>est -= X_bits.reshape(<NUM_LIT:1>, Ks.size, n_X, <NUM_LIT:1>)<EOL>est -= Y_bits.reshape(<NUM_LIT:1>, Ks.size, <NUM_LIT:1>, n_Y)<EOL>est /= <NUM_LIT:2><EOL>est += np.log(-<NUM_LIT:1> + X_ns[None, None, :, None] + Y_ns[None, None, None, :])<EOL>est += psi(Ks)[None, :, None, None]<EOL>if to_self:<EOL><INDENT>est[:, :, xrange(n_X), xrange(n_Y)] = <NUM_LIT:0><EOL><DEDENT>if clamp:  <EOL><INDENT>np.maximum(<NUM_LIT:0>, est, out=est)<EOL>np.minimum(np.log(<NUM_LIT:2>), est, out=est)<EOL><DEDENT>return est[:, :, :, :, None]<EOL>
r'''
    Estimate the difference between the Shannon entropy of an equally-weighted
    mixture between X and Y and the mixture of the Shannon entropies:

        JS(X, Y) = H[ (X + Y) / 2 ] - (H[X] + H[Y]) / 2

    We use a special case of the Hino-Murata weighted information estimator with
    a fixed M = n \alpha, about equivalent to the K-nearest-neighbor approach
    used for the other estimators:

        Hideitsu Hino and Noboru Murata (2013).
        Information estimators for weighted observations. Neural Networks.
        http://linkinghub.elsevier.com/retrieve/pii/S0893608013001676


    The estimator for JS(X, Y) is:

        log volume of the unit ball - log M + log(n + m - 1) + digamma(M)
        + 1/2 mean_X( d * log radius of largest ball in X+Y around X_i
                                with no more than M/(n+m-1) weight
                                where X points have weight 1 / (2 n - 1)
                                  and Y points have weight n / (m (2 n - 1))
                      - digamma(# of neighbors in that ball) )
        + 1/2 mean_Y( d * log radius of largest ball in X+Y around Y_i
                                with no more than M/(n+m-1) weight
                                where X points have weight m / (n (2 m - 1))
                                  and Y points have weight 1 / (2 m - 1)
                      - digamma(# of neighbors in that ball) )

        - 1/2 (log volume of the unit ball - log M + log(n - 1) + digamma(M))
        - 1/2 mean_X( d * log radius of the largest ball in X around X_i
                                with no more than M/(n-1) weight
                                where X points have weight 1 / (n - 1))
                      - digamma(# of neighbors in that ball) )

        - 1/2 (log volume of the unit ball - log M + log(m - 1) + digamma(M))
        - 1/2 mean_Y( d * log radius of the largest ball in Y around Y_i
                                with no more than M/(n-1) weight
                                where X points have weight 1 / (m - 1))
                      - digamma(# of neighbors in that ball) )

        =

        log(n + m - 1) + digamma(M)
        + 1/2 mean_X( d * log radius of largest ball in X+Y around X_i
                                with no more than M/(n+m-1) weight
                                where X points have weight 1 / (2 n - 1)
                                  and Y points have weight n / (m (2 n - 1))
                      - digamma(# of neighbors in that ball) )
        + 1/2 mean_Y( d * log radius of largest ball in X+Y around Y_i
                                with no more than M/(n+m-1) weight
                                where X points have weight m / (n (2 m - 1))
                                  and Y points have weight 1 / (2 m - 1)
                      - digamma(# of neighbors in that ball) )
        - 1/2 [log(n-1) + mean_X( d * log rho_M(X_i) )]
        - 1/2 [log(m-1) + mean_Y( d * log rho_M(Y_i) )]