def cw_sampling(X, y=None):

    def phi_sampling(s, D):
        return tf.pow((1.0 + ((4.0 * s) / ((2.0 * D) - 3))), (- 0.5))
    D = tf.cast(tf.shape(X)[1], tf.float32)
    N = tf.cast(tf.shape(X)[0], tf.float32)
    D_int = tf.cast(D, tf.int32)
    N_int = tf.cast(N, tf.int32)
    if (y is None):
        y = silverman_rule_of_thumb(N)
    YDistr = tf.contrib.distributions.MultivariateNormalDiag(loc=tf.zeros(D_int, tf.float32), scale_diag=tf.ones(D_int, tf.float32))
    Y = YDistr.sample(N_int)
    T = (1.0 / ((2.0 * N) * tf.sqrt((m.pi * y))))
    A0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(X, 1)), axis=2)
    A = tf.reduce_sum(phi_sampling((A0 / (4 * y)), D))
    B0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(Y, 0), tf.expand_dims(Y, 1)), axis=2)
    B = tf.reduce_sum(phi_sampling((B0 / (4 * y)), D))
    C0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(Y, 1)), axis=2)
    C = tf.reduce_sum(phi_sampling((C0 / (4 * y)), D))
    return (T * ((A + B) - (2 * C)))
def cw_sampling(X, y=None):

    def phi_sampling(s, D):
        return tf.pow((1.0 + ((4.0 * s) / ((2.0 * D) - 3))), (- 0.5))
    D = tf.cast(tf.shape(X)[1], tf.float32)
    N = tf.cast(tf.shape(X)[0], tf.float32)
    D_int = tf.cast(D, tf.int32)
    N_int = tf.cast(N, tf.int32)
    if (y is None):
        y = silverman_rule_of_thumb(N)
    YDistr = tf.contrib.distributions.MultivariateNormalDiag(loc=tf.zeros(D_int, tf.float32), scale_diag=tf.ones(D_int, tf.float32))
    Y = YDistr.sample(N_int)
    T = (1.0 / ((2.0 * N) * tf.sqrt((m.pi * y))))
    A0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(X, 1)), axis=2)
    A = tf.reduce_sum(phi_sampling((A0 / (4 * y)), D))
    B0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(Y, 0), tf.expand_dims(Y, 1)), axis=2)
    B = tf.reduce_sum(phi_sampling((B0 / (4 * y)), D))
    C0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(Y, 1)), axis=2)
    C = tf.reduce_sum(phi_sampling((C0 / (4 * y)), D))
    return (T * ((A + B) - (2 * C)))