# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # t-SNE

# t-SNE is a new award winning technique for dimension reduction and data visualization. t-SNE not only captures the local structure of the higher dimension but also preserves the global structures of the data like clusters. It has stunning ability to produce well-defined segregated clusters. t-SNE is based on stochastic neighbor embedding(SNE). t-SNE was developed to address some of the problems in SNE.

# ## SNE
# For understanding t-SNE, we need to know about stochastic neighbor embedding and its shortcomings. Stochastic neighbor embedding uses probablistic approach to embedd a high dimension dataset into lower dimension by preserving the nieghborhood structure of the dataset. A gaussian probability distribution centered on each point is defined over all the potential neighbors of this point. SNE aims to minimize the difference in probablity distribution in higher dimension and lower dimension.  
# For each object, $i$ and it's neighbor $j$, we compute a $P_{i|j}$ which reflects the probability that $j$ is neighbor of $i$    
# $\hspace{7em} P_{i|j} = \frac{\exp(-d_{ij}^2)}{\Sigma_{k\ne i}\exp(-d_{ij}^2)})$ where  
# $\hspace{7em} d_{ij}^2$ is the dissimilarity between element $i$ and $j$ given as input or calculated from the dataset provided.  
# The dissimilarity between $x_i$ and $x_j$ can be calculated using the following formula  
# $\hspace{7em} d_{ij}^2 = \frac{||x_i-x_j||^2}{2\sigma_i^2}$, where $\sigma_i$ generally calculated through a binary search by equating the entropy of the distribution centered at $x_i$ to perplexity which is chosen by hand. This method generates a probability matrix which is asymmetric.
#
#

# Next, SNE does the same calculation for every $Y_i$ and $Y_j$ in the lower dimension with gaussian probability distribution with $\sigma = 0.5$  
# $\hspace{7em} q_{ij} = \frac{\exp(-||y_i-y_j||^2)}{\Sigma_k \exp(-||y_k-y_i||^2)}$  
# Now, SNE tries to minimize the difference between these two distributions. We can calculate the difference between two distributions using Kullback-Liebler divergence. For two discrete distirbution $P$ and $Q$ KL divergence is given by  
# $\hspace{7em} D_{KL}(P||Q) = \Sigma_i P_{i}\frac{P_{i}}{Q_{i}}$.   
# SNE defines a cost function based of the difference between $p_{ij}$ and $q_{ij}$ which is given by  
# $\hspace{7em} C = \Sigma_i\Sigma_j P_{ij}log(\frac{P_{ij}}{q_{ij}})$  
# While embedding the dataset in lower dimension, two kinds of error can occur, first neighbors are mapped as faraway points($p_{ij}$ is large and  $q_{ij}$ is small) and points which are far away mapped as neighbors( $p_{ij}$ is small while $q_{ij}$ is large). Look closely at the cost function, the cost of the first kind of error i.e. mapping large $P_{ij}$ with small $q_{ij}$ is smaller than the cost while mapping small $p_{ij}$ as  large $q_{ij}$. The gradient to optimize the cost function is given by  
# $\hspace{7em} \frac{\delta C}{\delta Y} = 2\Sigma_i (P_{ij}-q_{ij}+P_{ji}-q_{ji})(y_i - y_j)$
#

# Some of the shorcomings of SNE approach are asymmetric probability matrix $P$, crowding problem. As pointed out earlier the probability matrix $P$ is asymmetric. Suppose a point $X_i$ is far away from other points, it's $P_{ij}$ will be very small for all $j$. So, It will have little effect on the cost function and embedding it correctly in the lower dimension will be hard. 
#

# Any n-dimensional euclidean space can have an object with n+1 or less equidistant vertices not more than that. Now, when the intrinsic dimension of a dataset is high say 20, and we are reducing it's dimensions from 100 to 2 or 3 our solution will be affected by crowding problem. The amount of space availble to map close points in 10 or 15 dimension will always be greater than the space availble in 2 or 3 dimensions. In order to map close points properly, moderately distant points will be pushed too far. This will eat the gaps in original clusters and it will look like single giant cluster. 

# We need to brush up few more topics before we move to t-SNE.  
# **Student-t distribution** --  Student-t distribution is a continuous symmetric probability distribution function with heavy tails. It is has only one parameter *degree of freedom*. As the *degree of freedom* increases it approaches the normal distribution function.  When *degree of freedom* =1, it takes the form of cauchy distribution function and its probability density function is given by  
# $\hspace{7em} f(t) = \frac{1}{\pi (1+t^2)}$  
# **Entropy** -- Entrophy is measure of the average information contained in a data. For a variable $x$ with $pdf$  $p(x)$  , it is given by  
# $\hspace{7em} H(x) = - \Sigma_{i}\,(\,p\,(x_i) \times  log_2(\,p(\,x_i\,)))$  
# **Perpexility** -- In information theory , perplexity measures how good a probability distribution predicts a sample. A low perplexity indicates that distribution function is good at predicting sample. It is given by  
# $\hspace{7em} Perpx(x) = 2^{H(x)}$, where $H(x)$ is the entropy of the distribution.  
#

# ## t-SNE algorithm.
#

# t-SNE differs from SNE in two ways, first it uses a student-t distribution to measure the similarity between points $Y_i$ and $Y_j$ in the lower dimension,secondly for the higher dimension it uses symmetric probability distribution such that $P_{ij} = P_{ji}$. Let's discuss the symmetric probability distribution first.  
# t-SNE defines the probability $P_{ij}$ as  
# $\hspace{7em} P_{ij} = \frac{P_{ij}+P_{ij}}{2n}$  
# This formulation makes sure that $\Sigma_jP_{ij} \gt \frac{1}{2n}$ for every $x_i$ and $x_i$ makes a significant contribution to the cost function. The gradient defined above becomes much simpler now  
# $\hspace{7em} \frac{\delta C}{\delta Y_i} = 4 \Sigma_{ij}(P_{ij}-q_{ij})(Y_i-Y_j)$   
# With student-t distribution having one degree one freedom , pairwise probability can be defined as  
# $\hspace{7em} q_{ij} = \frac{(1+||y_i-y_j||^2)^{-1}}{\Sigma_{k\ne i}(1+||y_i-y_k||^2)^{-1}}$  
# Now the gradient changes to   
# $\hspace{7em}\frac{\delta C}{\delta Y_i} = 4 \Sigma_{ij}(P_{ij}-q_{ij})(Y_i-Y_j)(1+||y_i-y_j||^2)^{-1}$ 

# ### psuedocode for t-SNE
# $\text{begin}$
# 1. compute pairwise similarity $P_{ij}$ for every $i$ and $j$.  
# 2. Update the probability with $P_{ij}$ = $\frac{P_{ij}+P_{ji}}{2n}$  
# 3. choose a random solution $Y_0 = (y_0,y_1,y_2,...y_n)$ where $y_i \in R^d$
# 4. While not done:  
#    $\hspace{1em}$compute pairwise similairites for $Y_0$  
#    $\hspace{1em}$compute the gradient $\frac {\delta C}{\delta y_i}$  
#    $\hspace{1em}$update the solution $y_i^t = y_i^{t-1} + \eta \frac{\delta C}{\delta y_i} + (y_i^{t-1} - y_i^{t-2})$  
#    $\hspace{1em}$if $t \gt \text{max_iter}$ :  
#    $\hspace{2em}$break  
#    $\hspace{1em}$else  
#    $\hspace{2em}$t = t+1  

# Let's implement the algorithm step by step
#

# #### Compute pariwise similarities
# For computing pairwise similarities we need to know the variance $\sigma_i$ for the gaussian centered at $x_i$. 
# One might think why not set a single value of $sigma_i$ for every $x_i$. The density of data is likely to vary, we need smaller $sigma_i$ for places with higher densities and bigger $\sigma_i$ for places where points are far away. The entropy of the gaussian distribution centered at $x_i$ increases as $\sigma_i$ increases. To get the $\sigma_i$ we need to perform a binary search such that perplexity of the gaussian distribution centered at $x_i$ is equal to the perplexity specified by the user.  
# Now, If you are thinking how perplexity fits into all this. You can think of perplexity as the smooth measure for the number of neighbors.  
# Then, compute the pairwise similarity as    
# $\hspace{7em} P_{i|j} = \frac{\exp(-d_{ij}^2)}{\Sigma_{k\ne i}\exp(-d_{ij}^2)})$  
# Once all the pairwise similarities have been calculated, update the similarity using the following rule  
# $\hspace{7em} P_{ii} = 0 $  
# $\hspace{7em} P_{ij} = \frac{P_{ij}+P_{ji}}{2n}$
#

# +
from sklearn import datasets

data = datasets.load_digits(n_class=6)
digits = data['data']
target = data['target']

# +
# define simiality function
from sklearn.metrics import pairwise as pw
import numpy as np

#Hbeta is the log of perplexity
#instead of calculating sigma , we are calculating beta = 1/(2*sigma^2)

def entropy(dist, beta):
    
    
    p = np.exp(-dist*beta)
    
    sum_p = p.sum()
    p = p/sum_p
    
    res = (p*np.log(p)).sum()
    
    return(-res, p)

def calculate_pi(dists, perplexity,i):
    
    beta_min = float("-inf")
    beta_max = float("inf")
    log_u = np.log(perplexity)
    dists = np.delete(dists,[i])
    func = lambda x: entropy(dists, x)
    
    beta = 1.0
    
    f1,p = func(beta)
    
    iter = 0
    diff = f1-log_u
    while(True):
        if np.abs(diff) < 1e-5:
            
            p = np.insert(p,i,0)
            return(p)
        
        if diff > 0.0:
        
            beta_min = beta
            if beta_max == np.inf or beta_max == -np.inf:
                beta *= 2.0
            else:
                beta = (beta + beta_max)/2.0
        else:
            beta_max = beta
            if beta_min == np.inf or beta_min == -np.inf:
                beta /= 2.0
    

            else:
                beta = (beta + beta_min) / 2.0
        
        iter += 1
        f1,p = func(beta)
        diff = f1-log_u
        if iter > 50:
            p = np.insert(p,i,0)
            
            return(p)

def get_Pmat(X, perplexity):
    n,p = X.shape
    P = np.zeros((n,n))
    D = pw.euclidean_distances(X, squared = True)
    
    for i in range(n):
    
        P[i] = calculate_pi(D[i], perplexity, i)

    # update the probabilities
    P = (P + P.T)
    
    P = P/P.sum()
    return(P)



# -

# ### calculation of low dimensional embedding
# Let's have a close look at the gradient obtained    
# $\hspace{7em}\frac{\delta C}{\delta y_i} = 4 \Sigma_{ij}(P_{ij}-q_{ij})(y_i-y_j)(1+||y_i-y_j||^2)^{-1}$   
# The t-SNE gradient is better in two ways, first, it induces heavy cost when dissimilar points are modelled by small pairwise distances in lower dimension and secondly, though the cost induced in first case is high it doesn't approach infinity due to which separate clusters are not pushed very far from each other.  
# Update equation for $y_i$ is given as  
# $\hspace{7em} y_i^t = y_i^{t-1} + \eta \frac{\delta C}{\delta y_i} + (y_i^{t-1} - y_i^{t-2})$  
# $\hspace{7em}$ where $t$ is the iteration.

# +
# generate a random sampled 

def q(D):
    
    D = 1/(1+D)
    
    np.fill_diagonal(D,0.0)
    
    q = D/D.sum()
    
    return(q)


def get_gradient(Y,P,q,D_low):
    
    n, d = Y.shape
    grad = np.zeros((n,d))
    
    D = 1/(1+D_low)
    temp = (P - q)*D
    
    for i in range(n):
        grad[i,:] = np.sum(temp[i,:].reshape(n,1)*(Y[i,:]-Y), axis=0)
    
    return(grad)


def cost_KL(P,q):
    
    n = P.shape[0]
    q[range(n),range(n)] = 1.00
    temp = P/q
    temp[temp==0.0]=1.0
    cost = (P*np.log(temp)).sum()
    return(cost)
    

np.random.seed(1000)

def t_sne(X,perplexity, n_dim):
    
    P = get_Pmat(X, perplexity)
    
    # value of lower dimension
    n,p = X.shape
    
    # generate a random sample of values for Y
    Y = np.random.rand(n,n_dim)

    # calculate the dissimilarity matrix for lower dimension points
    D_low = pw.euclidean_distances(Y, squared = True)
    minus = lambda x,y: (x-y).sum()
    alpha = 0.5
    max_iter=100

    Y_t_1 = np.zeros( (n, n_dim) )
    Y_t_2 = np.zeros( (n, n_dim) )
    eta = 0.9
    cost_old = 0.0
    q_ = q(D_low)

    cost_old = cost_KL(P,q_)
    #print("init cost = ",cost_old)
    
    for i in range(max_iter):
    
        grad = get_gradient(Y,P,q_,D_low)
    
        diff =  Y_t_1 - Y_t_2 
        Ynew = Y + eta*grad + alpha*(diff)
        Y_t_2 = Y_t_1
        Y_t_1 = Y
        Y = Ynew
    
        D_low = pw.euclidean_distances(Y, squared = True)
    
        q_ = q(D_low)
        cost_new = cost_KL(P,q_)
        #print("iter = ",i, "  cost = ",cost_new)
        if abs(cost_old - cost_new) <1e-4:
            break
        else:
            cost_old = cost_new

    return(Y)


# -

# How to select perplexity for the t-sne?
# Perplexity is an important parameter of the t-sne algorithm. We will understand it's effects on the toy dataset.
#

# +
# %matplotlib inline
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

target = datasets.load_iris().target

def plotter(x):
    
    Y = t_sne(iris, x,2)
    setosa = Y[target==0]
    versicolor = Y[target==1]
    verginica = Y[target==2]
    plt.scatter(setosa[:,0], setosa[:,1], c="b",label="setosa")
    plt.scatter(versicolor[:,0], versicolor[:,1], c="g",label="versicolor")
    plt.scatter(verginica[:,0], verginica[:,1], c="r",label="verginica")
    plt.legend()
    plt.show()

interact(plotter, x= widgets.IntSlider(min=5, max=50, value=10, step=5))
# -

# If you change value of perplexity from the slider you will see that, it's culsters changes dramatically. For t_sne to be meaningful we have to choose right value of perplexity. Perplexity balances the local and global aspects of the dataset. Very high value will lead to the merging of clusters into a single big clusters and low will produce many close small clusters which will be meaningless.

# ### t-SNE on pyspark

# Now, let's implement this algorithm on pyspark.

# +
# necessary imports
from sklearn import datasets
from pyspark.sql import SQLContext as SQC
from pyspark.ml.linalg import Vectors, VectorUDT
from pyspark.mllib.linalg.distributed import IndexedRowMatrix , IndexedRow

from pyspark.sql.types import *
from pyspark.sql.functions import udf, struct, collect_list, lit,col
import math as m
import numpy as np


# read the dataset
iris = datasets.load_iris()
iris_rdd = sc.parallelize(iris.data).zipWithIndex().map(lambda x:(x[1], Vectors.dense(x[0])))
iris_df  = sqc.createDataFrame(iris_rdd, ["id", "features"])
n,p = iris.data.shape
n_dim = 2



def weight(x, avg):
    
    if int(x)==1:
        return(0.0)
    elif x>2.5*avg:
        return(x)
    else:
        return(0.0)

    
def sort_list(x):
    
    y = sorted(x, key = lambda x: x[0])
    y = [s[1] for s in y]
    return(y)


# calculate the euclidean distance
udf_dist = udf(lambda x, y:  float(x.squared_distance(y)), DoubleType())

udf_weight = udf(weight, DoubleType())

# sort the distances according to id
udf_sort = udf(sort_list, ArrayType(DoubleType()))

# create the function to calculate the eucidean distances vector for each point
def dist_df(X_df):
    
    X_2 = X_df.crossJoin(X_df ).toDF('x_id', 'x_feature', 'y_id', 'y_feature')
    X_dist = X_2.withColumn("sim", udf_dist(X_2.x_feature, X_2.y_feature))
    avg = X_dist.groupBy().mean("sim").collect()[0][0]

    X_dist = X_dist.drop("x_feature","y_feature")

    nameCol = struct([name for name in ["y_id", "sim"]]).alias("map")
    X_dist = X_dist.select("x_id", nameCol)
    X_dist  = X_dist.groupby("x_id").agg(collect_list("map"))
    Dist = X_dist.withColumn("sims", udf_sort("collect_list(map)")).drop("collect_list(map)")
    
    return(Dist.select(col("x_id").alias("id"),col("sims").alias("Dist")))

DF_dist = dist_df(iris_df)

# entropy and p calculation is same as defined above

def entropy(dist, beta):
    
    
    p = np.exp(-dist*beta)
    
    sum_p = p.sum()
    p = p/sum_p
    
    res = (p*np.log(p)).sum()
    
    return(-res, p)

def get_p(dists, perplexity,i):
    dists = np.array(dists)
    beta_min = float("-inf")
    beta_max = float("inf")
    log_u = np.log(30)
    dists = np.delete(dists,[i])
    func = lambda x: entropy(dists, x)
    
    beta = 1.0
    
    f1,p = func(beta)
    
    iter = 0
    diff = f1-log_u
    while(True):
        if np.abs(diff) < 1e-5:
            
            p = np.insert(p,i,0)
            
            return(p.tolist())
        
        if diff > 0.0:
        
            beta_min = beta
            if beta_max == np.inf or beta_max == -np.inf:
                beta *= 2.0
            else:
                beta = (beta + beta_max)/2.0
        else:
            beta_max = beta
            if beta_min == np.inf or beta_min == -np.inf:
                beta /= 2.0
    

            else:
                beta = (beta + beta_min) / 2.0
        
        iter += 1
        f1,p = func(beta)
        diff = f1-log_u
        if iter > 50:
            p = np.insert(p,i,0)
            
            return(p.tolist())

udf_pi = udf(get_p, ArrayType(DoubleType()))

DF = DF_dist.withColumn("P", udf_pi("Dist",lit(30),"id")).drop("Dist")

# change pij to symmetric probability distribution

P_bm = IndexedRowMatrix(DF.select("id","P").rdd.map(lambda x: IndexedRow(x.id, x.P))).toBlockMatrix(10,10)
P = P_bm.add(P_bm.transpose())

udf_sum = udf(lambda x: float(np.sum(x)), DoubleType())

DF = (P.toIndexedRowMatrix().rows.map(lambda x:(x.index, x.vector) ).toDF(["id","P"])
                           .withColumn("Psum",udf_sum("P")))

# normalize the P matrix
Psum  = DF.groupBy().sum("Psum").collect()[0]["sum(Psum)"]

udf_divide = udf(lambda x,Psum:(np.array(x)/Psum).tolist(), ArrayType(DoubleType()))
DF = DF.withColumn("P", udf_divide("P",lit(Psum))).drop("Psum")

# ransom sample for the intial Y
np.random.seed(100)

Y = np.random.rand(n,n_dim).tolist()
Y_rdd = sc.parallelize(Y).zipWithIndex()

Y_df = spark.createDataFrame(Y_rdd,["Y","_id"])

DF = DF.join(Y_df,Y_df["_id"]==DF["id"]).drop("_id")

def sort_foo(x):
    
    x = sorted(x, key = lambda y:y[0])
    dists = list(map(lambda y:y[1][0], x))
    q = list(map(lambda y:1/(1 + y**2), dists))
    qsum  = float(np.sum(q))
    subs = list(map(lambda y:y[1][1], x))
    return([qsum, q, subs])

# function to calculate the q_ij from Yi and Yj

def Q_sub(df_Y):
    
    df2 = (df_Y.select(col("id").alias("id1"),col("features").alias("features1")).
            crossJoin(df_Y.select(col("id").alias("id2"),col("features").alias("features2"))))

    udf_dist = udf(lambda x, y:  float(np.sum((np.array(x)-np.array(y))**2)), DoubleType())
    udf_sub  = udf(lambda x,y:(np.array(x)-np.array(y)).tolist(), ArrayType(DoubleType()))

    df_dist = (df2.withColumn("sim", udf_dist(df2.features1, df2.features2))
            .withColumn("sub",udf_sub("features1","features2")))

    df3 = df_dist.select("id1","id2",struct("sim","sub").alias("val"))
    df4 = df3.select("id1",struct("id2","val").alias("set"))

    st = StructType([StructField("qsum", DoubleType()),
                     StructField("q", ArrayType(DoubleType())),
                 StructField("sub",ArrayType(ArrayType(DoubleType())))])

    udf_sort = udf(sort_foo,st)

    df5 = df4.groupBy("id1").agg(collect_list("set").alias("set"))
    df6 = (df5.withColumn("set",udf_sort("set")).select(col("id1").alias("_id"),
                                                col("set.qsum").alias("Qsum"),
                                                col("set.q").alias("Q"),
                                                col("set.sub").alias("Zsub")))
        
    return(df6)


Y = DF.select("id", col("Y").alias("features"))
udf_foo = udf(lambda x:np.zeros(x).tolist(), ArrayType(DoubleType()))

DF1 = Q_sub(Y)
DF = DF.join(DF1, DF["id"]==DF1["_id"]).drop("_id")
Qsum = DF.groupBy().sum("Qsum").collect()[0]["sum(Qsum)"]
DF = DF.withColumn("Y_t_1",udf_foo(lit(n_dim)))
DF = DF.select("id","P","Q","Zsub","Y", "Y_t_1",col("Y_t_1").alias("Y_t_2"))

# cost function

def kld_error(Q, P, Qsum):
    
    Q,P = np.array(Q)/Qsum, np.array(P)
    P[P==0.0] = 1.0
    cost = float(np.sum(P*np.log(P/Q)))
    return(cost)

udf_kld = udf(kld_error, DoubleType())

DF = DF.withColumn("err",udf_kld("Q","P",lit(Qsum)))

DF.cache()

error_old = DF.groupBy().sum("err").collect()[0]["sum(err)"]
DF = DF.drop("err")
#print("error-old = ",error_old)

#function for one descent step 
def descent_step(q, P, Y_t, Y_t_1, Y_t_2,Z_sub, Q_sum,alpha, eta):
    
    q,P,Y_t,Y_t_1,Y_t_2,Z_sub = list(map(lambda x:np.array(x), [q, P, Y_t, Y_t_1, Y_t_2, Z_sub] ) )
    Q_ = q/Q_sum 
    temp = ((P - Q_)*q).reshape(P.shape[0],1)
    dY = (temp*Z_sub).sum(axis=0)
    Y = Y_t + eta*dY +alpha*(Y_t_1-Y_t_2)
    
    return(Y.tolist())
    
udf_step = udf(descent_step, ArrayType(DoubleType()))

alpha =0.2
eta = 0.9


max_iter =10
for i in range(max_iter):
    DF = (DF.withColumn("new_Y", udf_step("Q","P","Y","Y_t_1",
                                          "Y_t_2","Zsub",lit(Qsum),lit(alpha),lit(eta))))
    
    DF = (DF.select("id","P",col("new_Y").alias("Y"),
                   col("Y").alias("Y_t_1"),
                   col("Y_t_1").alias("Y_t_2")))
    
    DF1 = DF.select("id",col("Y").alias("features"))
    DF1 = Q_sub(DF1)
    
    DF = DF.join(DF1, DF["id"]==DF1["_id"])
    Qsum = DF.groupBy().sum("Qsum").collect()[0]["sum(Qsum)"]
    
    DF = DF.withColumn("err",udf_kld("Q","P",lit(Qsum)))
    DF.cache()
    error_new = DF.groupBy().sum("err").collect()[0]["sum(err)"]
    DF = DF.drop("err")
    print("i = ",i,"  error = ",error_new)
    if abs(error_old - error_new)<1e-3:
        break
    else:
        error_old = error_new
        

DF.cache()
DF.show()
# -

# We have implemented a simple version of t-sne algorithm. There is a lot of scope for improvement both in algorithm and in gradient descent. There is a much faster version t-sne known as `berney-hut t-SNE`. It uses tree based algorithm to speed up searching. 

# Visualization of **MNIST** dataset

# ### Advantages and Drawbacks
#
# t-SNE works well for non-linear datasets. It work much better than other non-linear algorithms. Problems arise when intrinsic dimensions are higher i.e. more than 2-3 dimensions. t-SNE has tendency to get stuck in local optima. 

StarcoderdataJupyterScriptsDedupFiltered