content: # Sebastian Raschka, 2015 (http://sebastianraschka.com)
  # Python Machine Learning - Code Examples
  #
  # Chapter 5 - Compressing Data via Dimensionality Reduction
  #
  # S. Raschka. Python Machine Learning. Packt Publishing Ltd., 2015.
  # GitHub Repo: https://github.com/rasbt/python-machine-learning-book
  #
  # License: MIT
  # https://github.com/rasbt/python-machine-learning-book/blob/master/LICENSE.txt
  
  
  import pandas as pd
  import numpy as np
  from sklearn.preprocessing import StandardScaler
  from sklearn.decomposition import PCA
  import matplotlib.pyplot as plt
  from matplotlib.colors import ListedColormap
  from sklearn.linear_model import LogisticRegression
  from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
  from sklearn.datasets import make_moons
  from sklearn.datasets import make_circles
  from sklearn.decomposition import KernelPCA
  from scipy.spatial.distance import pdist, squareform
  from scipy import exp
  from scipy.linalg import eigh
  from matplotlib.ticker import FormatStrFormatter
  
  # for sklearn 0.18's alternative syntax
  from distutils.version import LooseVersion as Version
  from sklearn import __version__ as sklearn_version
  if Version(sklearn_version) < '0.18':
      from sklearn.grid_search import train_test_split
      from sklearn.lda import LDA
  else:
      from sklearn.model_selection import train_test_split
      from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Unsupervised dimensionality reduction'
        ' via principal component analysis')
  print(50 * '-')
  
  df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
                        'machine-learning-databases/wine/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']
  
  print('Wine data excerpt:\n\n:', df_wine.head())
  
  
  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, random_state=0)
  
  sc = StandardScaler()
  X_train_std = sc.fit_transform(X_train)
  X_test_std = sc.transform(X_test)
  
  cov_mat = np.cov(X_train_std.T)
  eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
  
  print('\nEigenvalues \n%s' % eigen_vals)
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Total and explained variance')
  print(50 * '-')
  
  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 components')
  plt.legend(loc='best')
  # plt.tight_layout()
  # plt.savefig('./figures/pca1.png', dpi=300)
  plt.show()
  
  #############################################################################
  print(50 * '=')
  print('Section: Feature Transformation')
  print(50 * '-')
  
  # 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(reverse=True)
  
  w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
                 eigen_pairs[1][1][:, np.newaxis]))
  print('Matrix W:\n', 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('./figures/pca2.png', dpi=300)
  plt.show()
  
  print('Dot product:\n', X_train_std[0].dot(w))
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Principal component analysis in scikit-learn')
  print(50 * '-')
  
  pca = PCA()
  X_train_pca = pca.fit_transform(X_train_std)
  print('Variance explained ratio:\n', 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 class samples
      for idx, cl in enumerate(np.unique(y)):
          plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                      alpha=0.8, c=cmap(idx),
                      marker=markers[idx], label=cl)
  
  
  lr = LogisticRegression()
  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('./figures/pca3.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('./figures/pca4.png', dpi=300)
  plt.show()
  
  pca = PCA(n_components=None)
  X_train_pca = pca.fit_transform(X_train_std)
  print('Explaind variance ratio:\n', pca.explained_variance_ratio_)
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Supervised data compression via linear discriminant analysis'
        ' - Computing the scatter matrices')
  print(50 * '-')
  
  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]))
  
  
  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]))
  
  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]))
  
  
  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[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]))
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Selecting linear discriminants for the new feature subspace')
  print(50 * '-')
  
  eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
  
  # 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 decreasing 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('./figures/lda1.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)
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Projecting samples onto the new feature space')
  print(50 * '-')
  
  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] * (-1),
                  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('./figures/lda2.png', dpi=300)
  plt.show()
  
  
  #############################################################################
  print(50 * '=')
  print('Section: LDA via scikit-learn')
  print(50 * '-')
  
  lda = LDA(n_components=2)
  X_train_lda = lda.fit_transform(X_train_std, y_train)
  
  lr = LogisticRegression()
  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/lda3.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/lda4.png', dpi=300)
  plt.show()
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Implementing a kernel principal component analysis in Python')
  print(50 * '-')
  
  
  def rbf_kernel_pca(X, gamma, n_components):
      """
      RBF kernel PCA implementation.
  
      Parameters
      ------------
      X: {NumPy ndarray}, shape = [n_samples, 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_samples, 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
      # numpy.eigh returns them in sorted order
      eigvals, eigvecs = eigh(K)
  
      # Collect the top k eigenvectors (projected samples)
      X_pc = np.column_stack((eigvecs[:, -i]
                              for i in range(1, n_components + 1)))
  
      return X_pc
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Example 1: Separating half-moon shapes')
  print(50 * '-')
  
  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('./figures/half_moon_1.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('./figures/half_moon_2.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')
  ax[0].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
  ax[1].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
  
  # plt.tight_layout()
  # plt.savefig('./figures/half_moon_3.png', dpi=300)
  plt.show()
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Example 2: Separating concentric circles')
  print(50 * '-')
  
  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('./figures/circles_1.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('./figures/circles_2.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('./figures/circles_3.png', dpi=300)
  plt.show()
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Projecting new data points')
  print(50 * '-')
  
  
  def rbf_kernel_pca(X, gamma, n_components):
      """
      RBF kernel PCA implementation.
  
      Parameters
      ------------
      X: {NumPy ndarray}, shape = [n_samples, 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_samples, 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
      # numpy.eigh returns them in sorted order
      eigvals, eigvecs = eigh(K)
  
      # Collect the top k eigenvectors (projected samples)
      alphas = np.column_stack((eigvecs[:, -i]
                                for i in range(1, n_components + 1)))
  
      # Collect the corresponding eigenvalues
      lambdas = [eigvals[-i] for i in range(1, n_components + 1)]
  
      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]
  print('New data point x_new:', x_new)
  
  x_proj = alphas[25]  # original projection
  print('Original projection x_proj:', 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)
  print('Reprojection x_reproj:', 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.legend(scatterpoints=1)
  
  # plt.tight_layout()
  # plt.savefig('./figures/reproject.png', dpi=300)
  plt.show()
  
  
  #############################################################################
  print(50 * '=')
  print('Section: Kernel principal component analysis in scikit-learn')
  print(50 * '-')
  
  
  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('./figures/scikit_kpca.png', dpi=300)
  plt.show()
  