{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nSubspace PCA (PCA without mean centering)\n=========================================\n\nThis example plots an animated gif showing how we can perform principle component analysis (PCA) without mean centering and obtain the same eigen vector.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib as mpl\nimport matplotlib.animation as animation\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nfrom matplotlib.animation import FuncAnimation\nfrom scipy.constants import golden as g\n\n\ndef dataset_fixed_cov():\n    '''Generate 1 Gaussians samples with the same covariance matrix'''\n    n, dim = 300, 2\n    np.random.seed(0)\n    C = np.array([[0., -0.3], [0.6, .3]])\n    X = np.dot(np.random.randn(n, dim), C) + [2., 2.]\n    return X\n\ndef proj_variance(X, vec):\n    return np.var(X @ vec.T)\n\ndef normalize(vec):\n    return vec / np.linalg.norm(vec)\n\ndef is_normalized(vec):\n    return np.linalg.norm(vec) == 1.\n\ndef unit_vector_from_rad(rad):\n    return np.array([np.cos(rad), np.sin(rad)])\n\n# Generate dataset\nX = dataset_fixed_cov()\n# Nomralise X\nX = np.apply_along_axis(normalize, axis=1, arr=(X-X.mean(axis=0)))\n\n# Calculate the direction that maximises the variance\n# with eigen decomposition\neig_vals, eig_vecs = np.linalg.eig(np.cov((X).T))\ntarget_phi = [vec for val, vec in sorted(zip(eig_vals, eig_vecs.T), reverse=True)][0]\n\n# calculate the angle of the target phi\ntarget_rad = np.angle(target_phi[0]+target_phi[1]*1j)\n\n# Predefine the number of time steps\nN = 300\n# Create N steps to \"solve\" for target\nrads = np.random.normal(loc=0, scale=np.pi, size=N) * np.geomspace(1, 2**-16, num=N) + target_rad\n\nplt.style.use('seaborn-dark')\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(g*6, 3))\n\n# Plot the scatters that persists (isn't redrawnstart_deg) \nax1.scatter(*X.T, c='blue', label='Target dataset') # Dataset\nax1.scatter(*X.mean(axis=0), c='red', label='Mean') # Mean \nax1.scatter(*[0,0], c='black', label='Origin') # Origin\nax1.quiver(*[0,0], *target_phi, angles='xy',scale_units='xy', scale=1, linestyle='--', alpha=0.6)\n# and init the quiver.\nQ = ax1.quiver(*[0,0,0,0], angles='xy',scale_units='xy', scale=1)\nax1.set_xlim(-2,2)\nax1.set_ylim(-2,2)\n\nx_data, y_data = [], []\nvl = ax2.axvline(0, 0, 1, linestyle='--', color='black', alpha=0.6)\nhl = ax2.axhline(0, 0, 1, linestyle='--', color='black', alpha=0.6)\nln, = ax2.plot(x_data, y_data,  'r.', alpha=0.2)\nax2.set_xlim(target_rad-np.pi/2, target_rad+np.pi/2)\nax2.set_ylim(0, 1)\n\nplots = [ln, Q, vl, hl]\n\ndef update_quiver(num, Q, phi, var):\n    fig.suptitle(f'step {num}')\n    Q.set_UVC(*phi)\n    ax1.set_title(f'Eigenvector: x={phi[0]:0.2f}, y={phi[0]:0.2f}')\n    return Q\n\ndef update_scatter(num, ln, var, vl, hl):\n    global x_data\n    global y_data \n    x_data += [num]\n    y_data += [var]\n    ln.set_data(x_data, y_data)\n    vl.set_data(num, [0, 2])\n    hl.set_data([0, 2], var)\n    ax2.set_title(f'J = {var:0.4f}')\n    return ln, vl, hl\n\ndef update(num, ln, Q, vl, hl):\n    phi = unit_vector_from_rad(rads[num])\n    var = proj_variance(X, phi)\n    # ln, Q = lnQ\n    Q = update_quiver(num, Q, phi, var)\n    ln, vl, hl = update_scatter(rads[num], ln, var, vl, hl)\n    return [ln, Q, vl, hl],\n            \n\nani = FuncAnimation(fig, update, fargs=(plots), frames=range(1,N),\n    interval=20, blit=False)\n\nplt.show()\n# ani.save('../docs/_static/pca/subspace_pca.gif',  writer='imagemagick', fps=60)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}