<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <h1 align="center"> t-Stochastic Nearest Neighbour Embeddings</h1>
t-SNE is an embedding technique that is commonly used to map high-dimensional data to 2 or 3 dimensions so that it can be viewed in scatterplots. The main idea behind t-SNE is to preserve points with small pairwise distances in the original data to the dimension reduced data.
t-SNE minimizes the divergence between two probability distributions
Step1:  We use the sklearn.datasets.load_digits method to load the MNIST data.
Step2:  This dataset contains data for 1797 images. Each image is an 8*8 matrix stored as a flat-packed array.
Next we combine the data and target into a single dataframe.
Step3:  Next we find out How many images we have per label.
Step4:  Next we scale mnist_df so that every feature has zero mean and unit variance.
Pairwise Distances, P and $\sigma_i$s
Step8:  From the scaled data, we must calculate the $P_{ij}$s. To do this we first calculate the pairwise distances between each pair of rows in the input data. For efficiency's sake, we use the sklearn.metrics.pairwise_distances library function.
Next, we start with a given purplexity target and then calculate the individual sigmas based on that.
Step11:  Now we're going to set up TensorFlow for the KLD minimization problem.
Step12:  Let's compare that against what the sklearn implementation gives us
Step13:  Appendix
Step14:  First we calculate Q using the direct iterative algorithm which requires iterating over rows and columns of y this gives us a reference to test our vectorized implementation for correctness.
Step15:  To calculate Q in a vectorized way, we note that
$D[i, j] = (y[i] - y[j]) (y[i] - y[j])^T = norm(y[i])^2 + norm(y[j])^2 - 2 \times dot(a[i], a[j])$
For the entire 2D array y, we can generalize this to
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import seaborn as sns
Explanation: <h1 align="center"> t-Stochastic Nearest Neighbour Embeddings</h1>
t-SNE is an embedding technique that is commonly used to map high-dimensional data to 2 or 3 dimensions so that it can be viewed in scatterplots. The main idea behind t-SNE is to preserve points with small pairwise distances in the original data to the dimension reduced data.
t-SNE minimizes the divergence between two probability distributions:  
a distribution that measures pairwise similarities of the input objects.
a distribution that measures pairwise similarities of the corresponding low-dimensional points in the embedding. 
Assume we are given a data set of (high-dimensional) input objects $D = {x_1, x_2, ... x_n}$ and a distance function $d(x_i, x_j)$ computes a distance between a pair of objects. The t-SNE algorithm learn learns an s (typically 2 or 3) dimensional embedding in which each object is represented by a point $y_i \in \mathbb{R}^s$.
To this end, t-SNE defines joint probabilities $p_{ij}$ that measure the pairwise similarity between objects $x_i$ and $x_j$
by symmetrizing two conditional probabilities:
$$\begin{align}
p_{j \mid i} \
 &= 0,\space if \space i = j \
 &= \dfrac {exp(-d(x_i, x_j)^2 / 2 \sigma_i^2)} {\sum_{k \ne i} exp(-d(x_i, x_k)^2 / 2 \sigma_i^2)} otherwise
\end{align}$$
and 
$$\begin{align}
p_{ij} = \dfrac {p_{i \mid j} + p_{j \mid i)}} {2n}
\end{align}$$
In  the  above  equation,  the  bandwidth  of  the  Gaussian  kernels, $\sigma_i$ is  set  in  such  a  way that the  perplexity of the conditional distribution $P_{ij}$ equals a predefined perplexity $u$. As a result,  the  optimal  value  of
$\sigma_i$ varies  per  object:  in  regions  of  the  data  space  with  a higher  data  density, $\sigma_i$ tends  to  be  smaller  than  in  regions  of  the  data  space  with  lower density.
In the s-dimensional embedding, the similarities between two points $y_i$ and $y_j$, the low-dimensional models of  $x_i$ and $x_j$ are measured using a normalized heavy-tailed kernel. Specifically, the embedding similarity $q_{ij}$ between the two points between $y_i$ and $y_j$ is computed as a normalized Student-t kernel with a single degree of freedom.
$$\begin{align}
q_{ij} \
 &= 0,\space if \space i = j \
 &= \dfrac {{(1 + \lVert y_i - y_j \rVert ^2})^{-1}} {{\sum_{k \ne l} (1 + \lVert y_k - y_l \rVert ^2})^{-1}} otherwise
\end{align}$$
The locations of the embedding points $y_i$ are determined by minimizing the Kullback-Leibler divergence between the joint distributions P and Q:
$$C(P, Q) = KLD(P \mid\mid Q) = \sum_{i \ne j} {p_{ij}} log (\dfrac{p_{ij}} {q_{ij}})$$ 
Resources
<a href="https://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf">https://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf</a>
In the code example below, we load classic MNIST dataset and reduce the dimensions to 2.
Getting the Data
End of explanation
from sklearn.datasets import load_digits
digits_data = load_digits()
from IPython.display import display
display(dir(digits_data))
display(digits_data.data.shape)
display(digits_data.target.shape)
Explanation: We use the sklearn.datasets.load_digits method to load the MNIST data.
End of explanation
mnist_df = pd.DataFrame(index=digits_data.target, data=digits_data.data)
mnist_df.head()
Explanation: This dataset contains data for 1797 images. Each image is an 8*8 matrix stored as a flat-packed array.
Next we combine the data and target into a single dataframe.
End of explanation
image_counts = mnist_df.groupby(mnist_df.index)[0].count()
ax = image_counts.plot(kind='bar', title='Image count per label in data')
Explanation: Next we find out How many images we have per label.
End of explanation
from sklearn.preprocessing import scale
mnist_df_scaled = pd.DataFrame(index=mnist_df.index, 
                               columns=mnist_df.columns, 
                               data=scale(mnist_df))
mnist_df_scaled.head()
Explanation: Next we scale mnist_df so that every feature has zero mean and unit variance.
Pairwise Distances, P and $\sigma_i$s
End of explanation
MACHINE_PRECISION = np.finfo(float).eps
from sklearn.metrics import pairwise_distances
def optimal_sigma(dist_i, i, target_entropy, n_iter=100, entropy_diff=1E-7):
    
    For the pairwise distances between the i-th feature vector and every other feature vector in the original dataset, 
    execute a binary search for ``sigma`` such the entropy of the conditional probability distribution
    ${P_i}$ equals ``target_entropy`` at ``entropy_diff`` precision. Return the optimal sigma.
    
    Assume that the distances are not squared.
    
    Execute at most ``n_iter`` searches. Raise ``ValueError`` if we haven't found a decent enough
    value of ``sigma``.
    
    Note that dist_i.loc[i] is the distance of the i-th feature vector to itself, i.e. 0. 
    
    assert dist_i.loc[i] == 0
    
    # initial value of sigma
    sigma = 1.0
    
    # initial left and right boundaries for the binary search
    sigma_min, sigma_max = -np.inf, np.inf
    
    for _ in range(1, n_iter+1):
        
        # Evaluate the Gaussian kernel with current sigma
        r = dist_i.pow(2).div(2 * (sigma ** 2))
        s = np.exp(-r)
        
        # Recall that p(j|i) = 0 if i = j
        s.loc[i] = 0
        p = s / s.sum()
        
        # the np.maximum trick below avoids taking log of very small (< MACHINE_PRECISION) numbers 
        # and ending up with -inf
        entropy = - p.dropna().dot(np.log(np.maximum(p.dropna(), MACHINE_PRECISION)))
        
        if np.fabs(target_entropy - entropy) <= entropy_diff:
            break            
                
        if entropy > target_entropy:
            # new boundary is [sigma_min, sigma]
            sigma_max = sigma
            # if sigma_min is still open
            if not np.isfinite(sigma_min):
                sigma *= 0.5        
        else:
            # new boundary is [sigma, sigma_max]
            sigma_min = sigma
           
            # if sigma_max is still open
            if not np.isfinite(sigma_max):                
                sigma *= 2.0
       
        # If both the left and right boundaries are closed, new sigma 
        # is the midpoint of sigma_min and sigma_max
        if np.all(np.isfinite([sigma_min, sigma_max])):
            sigma = (sigma_min + sigma_max) / 2
        
    else:
        raise ValueError("Unable to find a sigma after [{}] iterations that matches target entropy: [{}]".format(
            n_iter, target_entropy))        
        
    return sigma
def calc_optimal_sigmas(df, target_purplexity):
    
    From the DataFrame of feature vectors, ``df``, calculate pairwise distances and then find the optimal values
    for the Gaussian kernels for each conditional probability distribution {P_i}
    
    target_entropy = np.log(target_purplexity)
    paired_dists = pd.DataFrame(data=pairwise_distances(df.values, metric='l2'))
    
    optimal_sigmas = paired_dists.apply(lambda row: optimal_sigma(row, row.name, target_entropy), axis=1)
    
    # p_joint = (p_cond + p_cond.T) / (2 * df.shape[0])
    # return p_joint
    return paired_dists, optimal_sigmas
def calc_p(df, target_purplexity=30):
    
    Calculate the joint distribution of P_{ij} for the original input vectors x. 
    Assume ``pairwise_dist`` are squared. 
    
    paired_dists, optimal_sigmas = calc_optimal_sigmas(df, target_purplexity)
        
    exps = np.exp(-paired_dists)
    p_cond = exps.div(2 * optimal_sigmas.pow(2), axis=1)
    p_cond.values[np.diag_indices_from(p_cond)] = 0
    
    p_cond = p_cond.div(p_cond.sum(axis=1), axis=0)
    n_points = p_cond.shape[0]
    p_joint = (p_cond + p_cond.T) / (2 * n_points)
    
    return p_joint
p_joint = calc_p(mnist_df_scaled)
Explanation: From the scaled data, we must calculate the $P_{ij}$s. To do this we first calculate the pairwise distances between each pair of rows in the input data. For efficiency's sake, we use the sklearn.metrics.pairwise_distances library function.
Next, we start with a given purplexity target and then calculate the individual sigmas based on that.
End of explanation
import tensorflow as tf
display(tf.__version__)
def pairwise_dist(tf_y):
    Calculate pairwise distances between each pair of vectors in tf_y.
    tf_norms = tf.square(tf.norm(tf_y, axis=1))
     
    tf_r1 = tf.expand_dims(tf_norms, axis=1)
    tf_r2 = tf.expand_dims(tf_norms, axis=0)
    
    tf_y_dot_yT = tf.matmul(tf_y, tf_y, transpose_b=True)
    
    tf_dot = tf.cast(tf_y_dot_yT, dtype=tf.float32)
    tf_r = tf_r1 + tf_r2
    tf_d1 = tf_r - 2 * tf_dot
    
    return tf_d1
def calc_q(tf_y):
    
    Calculate the joint distribution of two embeddings y_i and y_j in tensorflow.   
    Call from inside an active tensorflow session only.    
    
    tf_pdist = pairwise_dist(tf_y)
    
    tf_d = 1 / (1 + tf_pdist)
    tf_d = tf.matrix_set_diag(tf_d, tf.zeros(tf.shape(tf_d)[0])) 
    
    tf_q = tf.div(tf_d, tf.reduce_sum(tf_d))
    
    return tf_q
embedding_size = 2
n_points = p_joint.shape[0]
losses = []
n_iter = 1000
loss_epsilon = 1E-8
learning_rate = 0.2
current_graph = tf.Graph()
with current_graph.as_default():
    
    # Placeholder for the joint distribution P of feature vectors in original space
    # This is a constant w.r.t the KLD minimization
    tf_p_joint = tf.placeholder(dtype=tf.float32, name='p_joint', shape=[n_points, n_points])
    
    # Feature vectors in the embedding space - initialized by sampling from random distribution
    tf_y = tf.Variable(name='y', validate_shape=False,
                        dtype=tf.float32,                        
                        initial_value=tf.random_normal([n_points, embedding_size]))
    
    # One step for iterative KLD minimization
    
    # calculate joint distribution Q of embeddings
    tf_q_joint = calc_q(tf_y)
    
    # Both P and Q have zeros in the diagonals. Since we want to calculate log{P/Q},
    # We temporarily replace the 1s with 0s, so the log of the diagonals are zeros
    # and they don't contribute to the KLD value.
    p_diag_1 = tf.matrix_set_diag(tf_p_joint, tf.ones(n_points))
    q_diag_1 = tf.matrix_set_diag(tf_q_joint, tf.ones(n_points))    
    tf_log_p_by_q = tf.log(tf.div(p_diag_1, q_diag_1))
    
    kld = tf.reduce_sum(tf.multiply(tf_p_joint, tf_log_p_by_q))
    
    optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate, name='Adam')
    train_op = optimizer.minimize(kld, name='KLD_minimization')
        
    with tf.Session() as sess:
        # initialize tensorflow variables
        init = tf.global_variables_initializer()
        sess.run(init)
        
        feed_dict = {tf_p_joint: p_joint.astype(np.float32).values}
        
        # run the optimization step n_iter times, breaking out if two successive steps
        # produce an absolute change in the value of the loss function <= loss_epsilon
        for i in range(1, n_iter+1):            
            _, loss_val = sess.run([train_op, kld], feed_dict=feed_dict)            
            losses.append(loss_val)
            
            if i % 100 == 0:
                print("After iteration: {}, loss: {}".format(i, loss_val))
            
            if len(losses) >= 2:
                last_loss = losses[-2]
                loss_delta = np.abs(last_loss-loss_val)
                if loss_delta < loss_epsilon:
                    print("Exiting after %s iterations, loss_delta [{}] <= loss_epsilon [{}".format(
                        n_iter, loss_delta, loss_epsilon))
                    break
        
        y_embeddings = sess.run(tf_y, feed_dict=feed_dict)
pd.Series(losses).rolling(10).mean().plot()
embeddings_df = pd.DataFrame(index=mnist_df_scaled.index, data=y_embeddings)
plot_source = embeddings_df.reset_index().rename(columns={
    'index': 'label', 0: 'x', 1:'y'})
fg = sns.FacetGrid(data=plot_source, hue='label', size=10)
fg.map(plt.scatter, 'x', 'y').add_legend()
Explanation: Now we're going to set up TensorFlow for the KLD minimization problem.
End of explanation
from sklearn.manifold import TSNE
# Extract the embeddings and convert into a DataFrame
sk_embedded = TSNE(n_components=2).fit_transform(mnist_df_scaled.values)
sk_embedded = pd.DataFrame(index=mnist_df_scaled.index, data=sk_embedded)
# Display
sk_embedded = sk_embedded.reset_index().rename(columns={'index': 'label', 0: 'x', 1:'y'})
fg = sns.FacetGrid(data=sk_embedded, hue='label', size=10)
fg.map(plt.scatter, 'x', 'y').add_legend()
Explanation: Let's compare that against what the sklearn implementation gives us:
End of explanation
y = pd.DataFrame(index=range(3), columns=range(5), data=np.random.uniform(1, 5, size=[3, 5]))
y
Explanation: Appendix: Vectorized Calculation of $Q_{ij}$ in TensorFlow$
End of explanation
Q_simple = pd.DataFrame(index=y.index, columns=y.index, data=0.0)
for i in range(0, y.shape[0]):
    for j in range(0, i):
        assert i != j, (i, j)            
        md = y.loc[i, :].sub(y.loc[j, :])
        d = 1 + np.linalg.norm(md)**2
        Q_simple.loc[i, j] = 1 / d
        Q_simple.loc[j, i] = 1 / d        
Q_simple
Explanation: First we calculate Q using the direct iterative algorithm which requires iterating over rows and columns of y this gives us a reference to test our vectorized implementation for correctness.
End of explanation
norms = y.apply(np.linalg.norm, axis=1).values
r1 = np.atleast_2d(norms**2)
r2 = r1.T
d1 = r1 + r2
d2 = d1 - 2 * np.dot(y, y.T)
d2 += 1
d3 = 1 / d2
d3[np.diag_indices_from(d3)] = 0
Q_vectorized = pd.DataFrame(d3)
Q_vectorized
from pandas.util.testing import assert_frame_equal
assert_frame_equal(Q_simple, Q_vectorized, check_less_precise=True)
Explanation: To calculate Q in a vectorized way, we note that
$D[i, j] = (y[i] - y[j]) (y[i] - y[j])^T = norm(y[i])^2 + norm(y[j])^2 - 2 \times dot(a[i], a[j])$
For the entire 2D array y, we can generalize this to:
D = np.atleast_2d(r) + np.atleast_2d(r).T - 2 * np.dot(y, y.T)
where r is (vector) of norms of each vector in y.
End of explanation
<END_TASK>