text: # coding: utf-8
  
  
  import pandas as pd
  from sklearn.model_selection import train_test_split
  from sklearn.preprocessing import StandardScaler
  import numpy as np
  import matplotlib.pyplot as plt
  from sklearn.decomposition import PCA
  from matplotlib.colors import ListedColormap
  from sklearn.linear_model import LogisticRegression
  from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
  from scipy.spatial.distance import pdist, squareform
  from scipy.linalg import eigh
  from distutils.version import LooseVersion as Version
  from scipy import __version__ as scipy_version
      from numpy import exp
      from scipy import exp
  from sklearn.datasets import make_moons
  from sklearn.datasets import make_circles
  from sklearn.decomposition import KernelPCA
  
  # *Python Machine Learning 3rd Edition* by [Sebastian Raschka](https://sebastianraschka.com), Packt Publishing Ltd. 2019
  # 
  # Code Repository: https://github.com/rasbt/python-machine-learning-book-3rd-edition
  # 
  # Code License: [MIT License](https://github.com/rasbt/python-machine-learning-book-3rd-edition/blob/master/LICENSE.txt)
  
  # # Python Machine Learning - Code Examples
  
  # # Chapter 5 - Compressing Data via Dimensionality Reduction
  
  # Note that the optional watermark extension is a small IPython notebook plugin that I developed to make the code reproducible. You can just skip the following line(s).
  
  
  
  
  
  # *The use of `watermark` is optional. You can install this Jupyter extension via*  
  # 
  #     conda install watermark -c conda-forge  
  # 
  # or  
  # 
  #     pip install watermark   
  # 
  # *For more information, please see: https://github.com/rasbt/watermark.*
  
  
  # ### Overview
  
  # - [Unsupervised dimensionality reduction via principal component analysis 128](#Unsupervised-dimensionality-reduction-via-principal-component-analysis-128)
  #   - [The main steps behind principal component analysis](#The-main-steps-behind-principal-component-analysis)
  #   - [Extracting the principal components step-by-step](#Extracting-the-principal-components-step-by-step)
  #   - [Total and explained variance](#Total-and-explained-variance)
  #   - [Feature transformation](#Feature-transformation)
  #   - [Principal component analysis in scikit-learn](#Principal-component-analysis-in-scikit-learn)
  # - [Supervised data compression via linear discriminant analysis](#Supervised-data-compression-via-linear-discriminant-analysis)
  #   - [Principal component analysis versus linear discriminant analysis](#Principal-component-analysis-versus-linear-discriminant-analysis)
  #   - [The inner workings of linear discriminant analysis](#The-inner-workings-of-linear-discriminant-analysis)
  #   - [Computing the scatter matrices](#Computing-the-scatter-matrices)
  #   - [Selecting linear discriminants for the new feature subspace](#Selecting-linear-discriminants-for-the-new-feature-subspace)
  #   - [Projecting examples onto the new feature space](#Projecting-examples-onto-the-new-feature-space)
  #   - [LDA via scikit-learn](#LDA-via-scikit-learn)
  # - [Using kernel principal component analysis for nonlinear mappings](#Using-kernel-principal-component-analysis-for-nonlinear-mappings)
  #   - [Kernel functions and the kernel trick](#Kernel-functions-and-the-kernel-trick)
  #   - [Implementing a kernel principal component analysis in Python](#Implementing-a-kernel-principal-component-analysis-in-Python)
  #     - [Example 1 – separating half-moon shapes](#Example-1:-Separating-half-moon-shapes)
  #     - [Example 2 – separating concentric circles](#Example-2:-Separating-concentric-circles)
  #   - [Projecting new data points](#Projecting-new-data-points)
  #   - [Kernel principal component analysis in scikit-learn](#Kernel-principal-component-analysis-in-scikit-learn)
  # - [Summary](#Summary)
  
  
  
  
  
  
  # # Unsupervised dimensionality reduction via principal component analysis
  
  # ## The main steps behind principal component analysis
  
  
  
  
  
  # ## Extracting the principal components step-by-step
  
  
  
  
  df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
                        'machine-learning-databases/wine/wine.data',
                        header=None)
  
  # if the Wine dataset is temporarily unavailable from the
  # UCI machine learning repository, un-comment the following line
  # of code to load the dataset from a local path:
  
  # df_wine = pd.read_csv('wine.data', header=None)
  
  df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
                     'Alcalinity of ash', 'Magnesium', 'Total phenols',
                     'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
                     'Color intensity', 'Hue',
                     'OD280/OD315 of diluted wines', 'Proline']
  
  df_wine.head()
  
  
  
  # Splitting the data into 70% training and 30% test subsets.
  
  
  
  
  X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
  
  X_train, X_test, y_train, y_test =     train_test_split(X, y, test_size=0.3, 
                       stratify=y,
                       random_state=0)
  
  
  # Standardizing the data.
  
  
  
  
  sc = StandardScaler()
  X_train_std = sc.fit_transform(X_train)
  X_test_std = sc.transform(X_test)
  
  
  # ---
  # 
  # **Note**
  # 
  # Accidentally, I wrote `X_test_std = sc.fit_transform(X_test)` instead of `X_test_std = sc.transform(X_test)`. In this case, it wouldn't make a big difference since the mean and standard deviation of the test set should be (quite) similar to the training set. However, as remember from Chapter 3, the correct way is to re-use parameters from the training set if we are doing any kind of transformation -- the test set should basically stand for "new, unseen" data.
  # 
  # My initial typo reflects a common mistake is that some people are *not* re-using these parameters from the model training/building and standardize the new data "from scratch." Here's simple example to explain why this is a problem.
  # 
  # Let's assume we have a simple training set consisting of 3 examples with 1 feature (let's call this feature "length"):
  # 
  # - train_1: 10 cm -> class_2
  # - train_2: 20 cm -> class_2
  # - train_3: 30 cm -> class_1
  # 
  # mean: 20, std.: 8.2
  # 
  # After standardization, the transformed feature values are
  # 
  # - train_std_1: -1.21 -> class_2
  # - train_std_2: 0 -> class_2
  # - train_std_3: 1.21 -> class_1
  # 
  # Next, let's assume our model has learned to classify examples with a standardized length value < 0.6 as class_2 (class_1 otherwise). So far so good. Now, let's say we have 3 unlabeled data points that we want to classify:
  # 
  # - new_4: 5 cm -> class ?
  # - new_5: 6 cm -> class ?
  # - new_6: 7 cm -> class ?
  # 
  # If we look at the "unstandardized "length" values in our training datast, it is intuitive to say that all of these examples are likely belonging to class_2. However, if we standardize these by re-computing standard deviation and and mean you would get similar values as before in the training set and your classifier would (probably incorrectly) classify examples 4 and 5 as class 2.
  # 
  # - new_std_4: -1.21 -> class 2
  # - new_std_5: 0 -> class 2
  # - new_std_6: 1.21 -> class 1
  # 
  # However, if we use the parameters from your "training set standardization," we'd get the values:
  # 
  # - example5: -18.37 -> class 2
  # - example6: -17.15 -> class 2
  # - example7: -15.92 -> class 2
  # 
  # The values 5 cm, 6 cm, and 7 cm are much lower than anything we have seen in the training set previously. Thus, it only makes sense that the standardized features of the "new examples" are much lower than every standardized feature in the training set.
  # 
  # ---
  
  # Eigendecomposition of the covariance matrix.
  
  
  
  cov_mat = np.cov(X_train_std.T)
  eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
  
  print('\nEigenvalues \n%s' % eigen_vals)
  
  
  # **Note**: 
  # 
  # Above, I used the [`numpy.linalg.eig`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html) function to decompose the symmetric covariance matrix into its eigenvalues and eigenvectors.
  #     <pre>>>> eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)</pre>
  #     This is not really a "mistake," but probably suboptimal. It would be better to use [`numpy.linalg.eigh`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html) in such cases, which has been designed for [Hermetian matrices](https://en.wikipedia.org/wiki/Hermitian_matrix). The latter always returns real  eigenvalues; whereas the numerically less stable `np.linalg.eig` can decompose nonsymmetric square matrices, you may find that it returns complex eigenvalues in certain cases. (S.R.)
  # 
  
  
  # ## Total and explained variance
  
  
  
  tot = sum(eigen_vals)
  var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
  cum_var_exp = np.cumsum(var_exp)
  
  
  
  
  
  
  plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
          label='Individual explained variance')
  plt.step(range(1, 14), cum_var_exp, where='mid',
           label='Cumulative explained variance')
  plt.ylabel('Explained variance ratio')
  plt.xlabel('Principal component index')
  plt.legend(loc='best')
  plt.tight_layout()
  # plt.savefig('images/05_02.png', dpi=300)
  plt.show()
  
  
  
  # ## Feature transformation
  
  
  
  # Make a list of (eigenvalue, eigenvector) tuples
  eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
                 for i in range(len(eigen_vals))]
  
  # Sort the (eigenvalue, eigenvector) tuples from high to low
  eigen_pairs.sort(key=lambda k: k[0], reverse=True)
  
  
  
  
  w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
                 eigen_pairs[1][1][:, np.newaxis]))
  print('Matrix W:\n', w)
  
  
  # **Note**
  # Depending on which version of NumPy and LAPACK you are using, you may obtain the Matrix W with its signs flipped. Please note that this is not an issue: If $v$ is an eigenvector of a matrix $\Sigma$, we have
  # 
  # $$\Sigma v = \lambda v,$$
  # 
  # where $\lambda$ is our eigenvalue,
  # 
  # 
  # then $-v$ is also an eigenvector that has the same eigenvalue, since
  # $$\Sigma \cdot (-v) = -\Sigma v = -\lambda v = \lambda \cdot (-v).$$
  
  
  
  X_train_std[0].dot(w)
  
  
  
  
  X_train_pca = X_train_std.dot(w)
  colors = ['r', 'b', 'g']
  markers = ['s', 'x', 'o']
  
  for l, c, m in zip(np.unique(y_train), colors, markers):
      plt.scatter(X_train_pca[y_train == l, 0], 
                  X_train_pca[y_train == l, 1], 
                  c=c, label=l, marker=m)
  
  plt.xlabel('PC 1')
  plt.ylabel('PC 2')
  plt.legend(loc='lower left')
  plt.tight_layout()
  # plt.savefig('images/05_03.png', dpi=300)
  plt.show()
  
  
  
  # ## Principal component analysis in scikit-learn
  
  # **NOTE**
  # 
  # The following four code cells has been added in addition to the content to the book, to illustrate how to replicate the results from our own PCA implementation in scikit-learn:
  
  
  
  
  pca = PCA()
  X_train_pca = pca.fit_transform(X_train_std)
  pca.explained_variance_ratio_
  
  
  
  
  plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')
  plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')
  plt.ylabel('Explained variance ratio')
  plt.xlabel('Principal components')
  
  plt.show()
  
  
  
  
  pca = PCA(n_components=2)
  X_train_pca = pca.fit_transform(X_train_std)
  X_test_pca = pca.transform(X_test_std)
  
  
  
  
  plt.scatter(X_train_pca[:, 0], X_train_pca[:, 1])
  plt.xlabel('PC 1')
  plt.ylabel('PC 2')
  plt.show()
  
  
  
  
  
  def plot_decision_regions(X, y, classifier, resolution=0.02):
  
      # setup marker generator and color map
      markers = ('s', 'x', 'o', '^', 'v')
      colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
      cmap = ListedColormap(colors[:len(np.unique(y))])
  
      # plot the decision surface
      x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
      x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
      xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                             np.arange(x2_min, x2_max, resolution))
      Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
      Z = Z.reshape(xx1.shape)
      plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
      plt.xlim(xx1.min(), xx1.max())
      plt.ylim(xx2.min(), xx2.max())
  
      # plot examples by class
      for idx, cl in enumerate(np.unique(y)):
          plt.scatter(x=X[y == cl, 0], 
                      y=X[y == cl, 1],
                      alpha=0.6, 
                      color=cmap(idx),
                      edgecolor='black',
                      marker=markers[idx], 
                      label=cl)
  
  
  # Training logistic regression classifier using the first 2 principal components.
  
  
  
  
  pca = PCA(n_components=2)
  X_train_pca = pca.fit_transform(X_train_std)
  X_test_pca = pca.transform(X_test_std)
  
  lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
  lr = lr.fit(X_train_pca, y_train)
  
  
  
  
  plot_decision_regions(X_train_pca, y_train, classifier=lr)
  plt.xlabel('PC 1')
  plt.ylabel('PC 2')
  plt.legend(loc='lower left')
  plt.tight_layout()
  # plt.savefig('images/05_04.png', dpi=300)
  plt.show()
  
  
  
  
  plot_decision_regions(X_test_pca, y_test, classifier=lr)
  plt.xlabel('PC 1')
  plt.ylabel('PC 2')
  plt.legend(loc='lower left')
  plt.tight_layout()
  # plt.savefig('images/05_05.png', dpi=300)
  plt.show()
  
  
  
  
  pca = PCA(n_components=None)
  X_train_pca = pca.fit_transform(X_train_std)
  pca.explained_variance_ratio_
  
  
  
  # # Supervised data compression via linear discriminant analysis
  
  # ## Principal component analysis versus linear discriminant analysis
  
  
  
  
  
  # ## The inner workings of linear discriminant analysis
  
  
  # ## Computing the scatter matrices
  
  # Calculate the mean vectors for each class:
  
  
  
  np.set_printoptions(precision=4)
  
  mean_vecs = []
  for label in range(1, 4):
      mean_vecs.append(np.mean(X_train_std[y_train == label], axis=0))
      print('MV %s: %s\n' % (label, mean_vecs[label - 1]))
  
  
  # Compute the within-class scatter matrix:
  
  
  
  d = 13 # number of features
  S_W = np.zeros((d, d))
  for label, mv in zip(range(1, 4), mean_vecs):
      class_scatter = np.zeros((d, d))  # scatter matrix for each class
      for row in X_train_std[y_train == label]:
          row, mv = row.reshape(d, 1), mv.reshape(d, 1)  # make column vectors
          class_scatter += (row - mv).dot((row - mv).T)
      S_W += class_scatter                          # sum class scatter matrices
  
  print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))
  
  
  # Better: covariance matrix since classes are not equally distributed:
  
  
  
  print('Class label distribution: %s' 
        % np.bincount(y_train)[1:])
  
  
  
  
  d = 13  # number of features
  S_W = np.zeros((d, d))
  for label, mv in zip(range(1, 4), mean_vecs):
      class_scatter = np.cov(X_train_std[y_train == label].T)
      S_W += class_scatter
  print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0],
                                                       S_W.shape[1]))
  
  
  # Compute the between-class scatter matrix:
  
  
  
  mean_overall = np.mean(X_train_std, axis=0)
  d = 13  # number of features
  S_B = np.zeros((d, d))
  for i, mean_vec in enumerate(mean_vecs):
      n = X_train_std[y_train == i + 1, :].shape[0]
      mean_vec = mean_vec.reshape(d, 1)  # make column vector
      mean_overall = mean_overall.reshape(d, 1)  # make column vector
      S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
  
  print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))
  
  
  
  # ## Selecting linear discriminants for the new feature subspace
  
  # Solve the generalized eigenvalue problem for the matrix $S_W^{-1}S_B$:
  
  
  
  eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
  
  
  # **Note**:
  #     
  # Above, I used the [`numpy.linalg.eig`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html) function to decompose the symmetric covariance matrix into its eigenvalues and eigenvectors.
  #     <pre>>>> eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)</pre>
  #     This is not really a "mistake," but probably suboptimal. It would be better to use [`numpy.linalg.eigh`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html) in such cases, which has been designed for [Hermetian matrices](https://en.wikipedia.org/wiki/Hermitian_matrix). The latter always returns real  eigenvalues; whereas the numerically less stable `np.linalg.eig` can decompose nonsymmetric square matrices, you may find that it returns complex eigenvalues in certain cases. (S.R.)
  # 
  
  # Sort eigenvectors in descending order of the eigenvalues:
  
  
  
  # Make a list of (eigenvalue, eigenvector) tuples
  eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
                 for i in range(len(eigen_vals))]
  
  # Sort the (eigenvalue, eigenvector) tuples from high to low
  eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
  
  # Visually confirm that the list is correctly sorted by decreasing eigenvalues
  
  print('Eigenvalues in descending order:\n')
  for eigen_val in eigen_pairs:
      print(eigen_val[0])
  
  
  
  
  tot = sum(eigen_vals.real)
  discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
  cum_discr = np.cumsum(discr)
  
  plt.bar(range(1, 14), discr, alpha=0.5, align='center',
          label='Individual "discriminability"')
  plt.step(range(1, 14), cum_discr, where='mid',
           label='Cumulative "discriminability"')
  plt.ylabel('"Discriminability" ratio')
  plt.xlabel('Linear discriminants')
  plt.ylim([-0.1, 1.1])
  plt.legend(loc='best')
  plt.tight_layout()
  # plt.savefig('images/05_07.png', dpi=300)
  plt.show()
  
  
  
  
  w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
                eigen_pairs[1][1][:, np.newaxis].real))
  print('Matrix W:\n', w)
  
  
  
  # ## Projecting examples onto the new feature space
  
  
  
  X_train_lda = X_train_std.dot(w)
  colors = ['r', 'b', 'g']
  markers = ['s', 'x', 'o']
  
  for l, c, m in zip(np.unique(y_train), colors, markers):
      plt.scatter(X_train_lda[y_train == l, 0],
                  X_train_lda[y_train == l, 1] * (-1),
                  c=c, label=l, marker=m)
  
  plt.xlabel('LD 1')
  plt.ylabel('LD 2')
  plt.legend(loc='lower right')
  plt.tight_layout()
  # plt.savefig('images/05_08.png', dpi=300)
  plt.show()
  
  
  
  # ## LDA via scikit-learn
  
  
  
  
  lda = LDA(n_components=2)
  X_train_lda = lda.fit_transform(X_train_std, y_train)
  
  
  
  
  
  lr = LogisticRegression(multi_class='ovr', random_state=1, solver='lbfgs')
  lr = lr.fit(X_train_lda, y_train)
  
  plot_decision_regions(X_train_lda, y_train, classifier=lr)
  plt.xlabel('LD 1')
  plt.ylabel('LD 2')
  plt.legend(loc='lower left')
  plt.tight_layout()
  # plt.savefig('images/05_09.png', dpi=300)
  plt.show()
  
  
  
  
  X_test_lda = lda.transform(X_test_std)
  
  plot_decision_regions(X_test_lda, y_test, classifier=lr)
  plt.xlabel('LD 1')
  plt.ylabel('LD 2')
  plt.legend(loc='lower left')
  plt.tight_layout()
  # plt.savefig('images/05_10.png', dpi=300)
  plt.show()
  
  
  
  # # Using kernel principal component analysis for nonlinear mappings
  
  
  
  
  
  
  # ## Implementing a kernel principal component analysis in Python
  
  
  
  
  
  
  
  if scipy_version >= Version('1.4.1'):
  else:
  
  
  
  
  def rbf_kernel_pca(X, gamma, n_components):
      """
      RBF kernel PCA implementation.
  
      Parameters
      ------------
      X: {NumPy ndarray}, shape = [n_examples, n_features]
          
      gamma: float
        Tuning parameter of the RBF kernel
          
      n_components: int
        Number of principal components to return
  
      Returns
      ------------
       X_pc: {NumPy ndarray}, shape = [n_examples, k_features]
         Projected dataset   
  
      """
      # Calculate pairwise squared Euclidean distances
      # in the MxN dimensional dataset.
      sq_dists = pdist(X, 'sqeuclidean')
  
      # Convert pairwise distances into a square matrix.
      mat_sq_dists = squareform(sq_dists)
  
      # Compute the symmetric kernel matrix.
      K = exp(-gamma * mat_sq_dists)
  
      # Center the kernel matrix.
      N = K.shape[0]
      one_n = np.ones((N, N)) / N
      K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
  
      # Obtaining eigenpairs from the centered kernel matrix
      # scipy.linalg.eigh returns them in ascending order
      eigvals, eigvecs = eigh(K)
      eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]
  
      # Collect the top k eigenvectors (projected examples)
      X_pc = np.column_stack([eigvecs[:, i]
                              for i in range(n_components)])
  
      return X_pc
  
  
  
  # ### Example 1: Separating half-moon shapes
  
  
  
  
  X, y = make_moons(n_samples=100, random_state=123)
  
  plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='^', alpha=0.5)
  plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='o', alpha=0.5)
  
  plt.tight_layout()
  # plt.savefig('images/05_12.png', dpi=300)
  plt.show()
  
  
  
  
  
  scikit_pca = PCA(n_components=2)
  X_spca = scikit_pca.fit_transform(X)
  
  fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
  
  ax[0].scatter(X_spca[y == 0, 0], X_spca[y == 0, 1],
                color='red', marker='^', alpha=0.5)
  ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
                color='blue', marker='o', alpha=0.5)
  
  ax[1].scatter(X_spca[y == 0, 0], np.zeros((50, 1)) + 0.02,
                color='red', marker='^', alpha=0.5)
  ax[1].scatter(X_spca[y == 1, 0], np.zeros((50, 1)) - 0.02,
                color='blue', marker='o', alpha=0.5)
  
  ax[0].set_xlabel('PC1')
  ax[0].set_ylabel('PC2')
  ax[1].set_ylim([-1, 1])
  ax[1].set_yticks([])
  ax[1].set_xlabel('PC1')
  
  plt.tight_layout()
  # plt.savefig('images/05_13.png', dpi=300)
  plt.show()
  
  
  
  
  X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)
  
  fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
  ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], 
              color='red', marker='^', alpha=0.5)
  ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],
              color='blue', marker='o', alpha=0.5)
  
  ax[1].scatter(X_kpca[y==0, 0], np.zeros((50, 1))+0.02, 
              color='red', marker='^', alpha=0.5)
  ax[1].scatter(X_kpca[y==1, 0], np.zeros((50, 1))-0.02,
              color='blue', marker='o', alpha=0.5)
  
  ax[0].set_xlabel('PC1')
  ax[0].set_ylabel('PC2')
  ax[1].set_ylim([-1, 1])
  ax[1].set_yticks([])
  ax[1].set_xlabel('PC1')
  
  plt.tight_layout()
  # plt.savefig('images/05_14.png', dpi=300)
  plt.show()
  
  
  
  # ### Example 2: Separating concentric circles
  
  
  
  
  X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)
  
  plt.scatter(X[y == 0, 0], X[y == 0, 1], color='red', marker='^', alpha=0.5)
  plt.scatter(X[y == 1, 0], X[y == 1, 1], color='blue', marker='o', alpha=0.5)
  
  plt.tight_layout()
  # plt.savefig('images/05_15.png', dpi=300)
  plt.show()
  
  
  
  
  scikit_pca = PCA(n_components=2)
  X_spca = scikit_pca.fit_transform(X)
  
  fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
  
  ax[0].scatter(X_spca[y == 0, 0], X_spca[y == 0, 1],
                color='red', marker='^', alpha=0.5)
  ax[0].scatter(X_spca[y == 1, 0], X_spca[y == 1, 1],
                color='blue', marker='o', alpha=0.5)
  
  ax[1].scatter(X_spca[y == 0, 0], np.zeros((500, 1)) + 0.02,
                color='red', marker='^', alpha=0.5)
  ax[1].scatter(X_spca[y == 1, 0], np.zeros((500, 1)) - 0.02,
                color='blue', marker='o', alpha=0.5)
  
  ax[0].set_xlabel('PC1')
  ax[0].set_ylabel('PC2')
  ax[1].set_ylim([-1, 1])
  ax[1].set_yticks([])
  ax[1].set_xlabel('PC1')
  
  plt.tight_layout()
  # plt.savefig('images/05_16.png', dpi=300)
  plt.show()
  
  
  
  
  X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)
  
  fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))
  ax[0].scatter(X_kpca[y == 0, 0], X_kpca[y == 0, 1],
                color='red', marker='^', alpha=0.5)
  ax[0].scatter(X_kpca[y == 1, 0], X_kpca[y == 1, 1],
                color='blue', marker='o', alpha=0.5)
  
  ax[1].scatter(X_kpca[y == 0, 0], np.zeros((500, 1)) + 0.02,
                color='red', marker='^', alpha=0.5)
  ax[1].scatter(X_kpca[y == 1, 0], np.zeros((500, 1)) - 0.02,
                color='blue', marker='o', alpha=0.5)
  
  ax[0].set_xlabel('PC1')
  ax[0].set_ylabel('PC2')
  ax[1].set_ylim([-1, 1])
  ax[1].set_yticks([])
  ax[1].set_xlabel('PC1')
  
  plt.tight_layout()
  # plt.savefig('images/05_17.png', dpi=300)
  plt.show()
  
  
  
  # ## Projecting new data points
  
  
  
  
  def rbf_kernel_pca(X, gamma, n_components):
      """
      RBF kernel PCA implementation.
  
      Parameters
      ------------
      X: {NumPy ndarray}, shape = [n_examples, n_features]
          
      gamma: float
        Tuning parameter of the RBF kernel
          
      n_components: int
        Number of principal components to return
  
      Returns
      ------------
       alphas: {NumPy ndarray}, shape = [n_examples, k_features]
         Projected dataset 
       
       lambdas: list
         Eigenvalues
  
      """
      # Calculate pairwise squared Euclidean distances
      # in the MxN dimensional dataset.
      sq_dists = pdist(X, 'sqeuclidean')
  
      # Convert pairwise distances into a square matrix.
      mat_sq_dists = squareform(sq_dists)
  
      # Compute the symmetric kernel matrix.
      K = exp(-gamma * mat_sq_dists)
  
      # Center the kernel matrix.
      N = K.shape[0]
      one_n = np.ones((N, N)) / N
      K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
  
      # Obtaining eigenpairs from the centered kernel matrix
      # scipy.linalg.eigh returns them in ascending order
      eigvals, eigvecs = eigh(K)
      eigvals, eigvecs = eigvals[::-1], eigvecs[:, ::-1]
  
      # Collect the top k eigenvectors (projected examples)
      alphas = np.column_stack([eigvecs[:, i]
                                for i in range(n_components)])
  
      # Collect the corresponding eigenvalues
      lambdas = [eigvals[i] for i in range(n_components)]
  
      return alphas, lambdas
  
  
  
  
  X, y = make_moons(n_samples=100, random_state=123)
  alphas, lambdas = rbf_kernel_pca(X, gamma=15, n_components=1)
  
  
  
  
  x_new = X[25]
  x_new
  
  
  
  
  x_proj = alphas[25] # original projection
  x_proj
  
  
  
  
  def project_x(x_new, X, gamma, alphas, lambdas):
      pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
      k = np.exp(-gamma * pair_dist)
      return k.dot(alphas / lambdas)
  
  # projection of the "new" datapoint
  x_reproj = project_x(x_new, X, gamma=15, alphas=alphas, lambdas=lambdas)
  x_reproj 
  
  
  
  
  plt.scatter(alphas[y == 0, 0], np.zeros((50)),
              color='red', marker='^', alpha=0.5)
  plt.scatter(alphas[y == 1, 0], np.zeros((50)),
              color='blue', marker='o', alpha=0.5)
  plt.scatter(x_proj, 0, color='black',
              label='Original projection of point X[25]', marker='^', s=100)
  plt.scatter(x_reproj, 0, color='green',
              label='Remapped point X[25]', marker='x', s=500)
  plt.yticks([], [])
  plt.legend(scatterpoints=1)
  
  plt.tight_layout()
  # plt.savefig('images/05_18.png', dpi=300)
  plt.show()
  
  
  
  # ## Kernel principal component analysis in scikit-learn
  
  
  
  
  X, y = make_moons(n_samples=100, random_state=123)
  scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
  X_skernpca = scikit_kpca.fit_transform(X)
  
  plt.scatter(X_skernpca[y == 0, 0], X_skernpca[y == 0, 1],
              color='red', marker='^', alpha=0.5)
  plt.scatter(X_skernpca[y == 1, 0], X_skernpca[y == 1, 1],
              color='blue', marker='o', alpha=0.5)
  
  plt.xlabel('PC1')
  plt.ylabel('PC2')
  plt.tight_layout()
  # plt.savefig('images/05_19.png', dpi=300)
  plt.show()
  
  
  
  # # Summary
  
  # ...
  
  # ---
  # 
  # Readers may ignore the next cell.
  
  
  
  
  