code: def fit(self, X):
          
  
          # compute the std. deviation with mean vector x in X
          self._compute_std_dev(X)
  
          # Kullback–Leibler divergence
          kl_cost = KL_Divergence()
  
          # compute high-dimensional affinities (Gaussian Distribution)
          high_dim_dist = self._get_high_dim_dist(X)
          # sample initial solutions
          Y = np.random.randn(X.shape[0], self.n_components)
  
          prev_Ys = [Y, Y]
  
          for iteration in range(1, self.n_iter+1):
              # compute low-dimensional affinities (Student t-Distribution)
              low_dim_dist = self._get_low_dim_dist(Y)
      
              for i in range(Y.shape[0]):
                  # compute gradient
                  grad = kl_cost.gradient(high_dim_dist, low_dim_dist, Y, i)
                  # set new Y[i]
                  Y[i] = prev_Ys[1][i] + self.learning_rate * grad + self.momentum * (prev_Ys[1][i] - prev_Ys[0][i])
  
              prev_Ys = [prev_Ys[1], Y]
  
              if iteration % 100 == 0 and self.verbose:
                  low_dim_dist = self._get_low_dim_dist(Y)
                  print(f"ITERATION: {iteration}{3*' '}|||{3*' '}KL divergence: {kl_cost(high_dim_dist, low_dim_dist)}")
  
          self.embeddings = Y
          return self
docs: Gradient Descent optimization process
          
          Tunes the embeddings (Y) so that their pairwise distance distribution
          matches the input high-dimensional data (X) pairwise distance distribution.
          In other words, minimizes the KL divergence cost.