A4.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X_qFHr9IhitQ"
      },
      "source": [
        "# **Machine Learning Assignment 4: Dimensionality Reduction**\n",
        "\n",
        "Read this notebook and do the exercises. Then:\n",
        "\n",
        "* Copy your exercise solutions into the skeleton code in `A4.py`, and test them with `pytest` and `A4_test.py`. \n",
        "\n",
        "When you're satisfied that your code in `A4.py` passes the tests, commit and push it back to GitHub. You can also push your notebook -- that will never negatively impact your grade, but could give us an idea of what you did wrong for open questions.\n",
        "\n",
        "* Make a PDF report (using the provided template) with answers to the open questions (as well as your graphics). It should be named `A4_s1234567890.pdf`, but using your own student number instead of 1234567890. Hand the PDF report in through BrightSpace. Your answers don't need to be long, just be to the point - _illustrative examples are welcome._\n",
        "\n",
        "* If you do not submit both `A4.py` to GitHub and the PDF report to BrightSpace, you will fail this assignment.\n",
        "\n",
        "* This is an individual assignment. Your code and report must be your own work.\n",
        "\n",
        "---\n",
        "\n",
        "Specific notes for this assignment:\n",
        "\n",
        "* Make sure to look at the code that is not part of programming exercises, we sometimes give it to you because you need it later on.\n",
        "\n",
        "* There are 8 programming exercises that are labelled with subsubsections that start with \"Exercise [number]:\".\n",
        "\n",
        "* There are 5 regular open questions, 2 plot questions, and 1 open bonus question that you should answer in your report. These are labelled with subsubsections that start with \"Question [number] [(bonus)  or (2 points) if applicable]:\".\n",
        "\n",
        "* The 8 programming exercises (`A4.py`) make up 40% of your grade. The 8 or 9 questions that end up in your report (`A4_[your ID].pdf`) make up the other 60% of your grade."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EcbMOtMvwSnl"
      },
      "source": [
        "# Set-Up"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jwnovVS-baoG"
      },
      "source": [
        "## Import packages and functions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-DLYyhnkwSnw"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.preprocessing import scale as sk_scale\n",
        "from sklearn.datasets import load_breast_cancer\n",
        "from sklearn.svm import SVC\n",
        "from sklearn.utils import shuffle as sk_shuffle"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KdW8IqR7bP6w"
      },
      "source": [
        "## Define some basic functions and classes\n",
        "\n",
        "The auxiliary functions and the `NFolds` class will be used bellow."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "kYnzzLgNymPR"
      },
      "outputs": [],
      "source": [
        "def helpful_eq(a, b, failing_is_good=False):\n",
        "    \"\"\"Basically `==` after rounding with prints. \n",
        "    `a` and `b` can be numbers or (>1D) (NumPy) arrays.\"\"\"\n",
        "    def print_bad_news():\n",
        "        print(a)\n",
        "        print(\"helpful_eq(...) fail: ^ does not equal:\")\n",
        "        print(b)\n",
        "\n",
        "    try:\n",
        "        if hasattr(a, \"__len__\"):\n",
        "            r = np.allclose(a, b, atol=1e-3)\n",
        "        else:\n",
        "            r = round(a, 3) == round(b, 3)\n",
        "        if failing_is_good:\n",
        "            r = not r\n",
        "        if not r:\n",
        "            print_bad_news()\n",
        "        return r\n",
        "    except Exception as e:\n",
        "        print_bad_news()\n",
        "        print(\"And/or we encountered exception message\")\n",
        "        print(e)\n",
        "        return False\n",
        "\n",
        "\n",
        "def mean_squared_error(true, pred):\n",
        "    \"\"\"`true` and `pred` should be (numpy.nd)arrays with similar shapes.\n",
        "    Returns mean (over instances/rows) of sum of squared feature difference (in columns).\"\"\"\n",
        "    return np.mean(((true - pred) ** 2).mean(axis=1))\n",
        "\n",
        "\n",
        "class NFolds:\n",
        "    def __init__(self, X, y, n_folds=5, seed=42):\n",
        "        \"\"\" Initialize the KFolds instance\n",
        "\n",
        "        :param X: numpy.ndarray of feature columns \n",
        "        :param y: numpy.ndarray of labels\n",
        "        :param n_folds: number of folds desired\n",
        "        :param seed: random seed, if you want reproducible results (optional)\n",
        "\n",
        "        After initialization, self.folds will store n_folds folds.\n",
        "        Each fold is a pair of arrays with training indices and test indices.\n",
        "        The folds are as evenly distributed in size as possible.\n",
        "        All the test segments are pairwise disjoint.\n",
        "        \"\"\"\n",
        "        self.X = X\n",
        "        self.y = y\n",
        "        self.n_folds = n_folds\n",
        "        self.folds = []\n",
        "        indices = np.arange(X.shape[0])\n",
        "        if seed is not None:\n",
        "            np.random.seed(seed=seed)\n",
        "        np.random.shuffle(indices)\n",
        "        fold_size = X.shape[0] / n_folds\n",
        "        for fold_num in range(n_folds):\n",
        "            test = indices[int(fold_num * fold_size): int((fold_num + 1) * fold_size)]\n",
        "            train = np.concatenate([indices[: int(fold_num * fold_size)],\n",
        "                                    indices[int((fold_num + 1) * fold_size):]])\n",
        "            self.folds.append((train, test))\n",
        "\n",
        "    def get_fold(self, fold_num):\n",
        "        \"\"\" Get the training and test data of the fold_num-th fold\n",
        "\n",
        "        :param fold_num: Which fold's division of the data to use\n",
        "        :return: Training and test features/labels\n",
        "        \"\"\"\n",
        "        train, test = self.folds[fold_num]\n",
        "        X_train = self.X[train]\n",
        "        X_test = self.X[test]\n",
        "        y_train = self.y[train]\n",
        "        y_test = self.y[test]\n",
        "        return X_train, X_test, y_train, y_test"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s_X3cHfqwSnx"
      },
      "source": [
        "# Principal component analysis\n",
        "## Matrix transformations and decompositions\n",
        "We can view feature construction and selection from a geometric perspective, assuming quantitative features. To this end we represent our data set as a matrix $\\mathbf{X}$ with $n$ data points in rows and $d$ features in columns, which we want to transform into a new matrix $\\mathbf{W}$ with $n$ rows and $k$ columns by means of matrix operations. To simplify matters a bit, we assume that $\\mathbf{X}$ is zero-centred and that $\\mathbf{W} = \\mathbf{XT}$ for some $d$-by-$k$ transformation matrix $\\mathbf{T}$. For example, feature scaling corresponds to $\\mathbf{T}$ being a $d$-by-$d$ diagonal matrix; this can be combined with feature selection by removing some of $\\mathbf{T}$’s columns. A rotation is achieved by $\\mathbf{T}$ being orthogonal, i.e., $\\mathbf{T} \\mathbf{T}^{\\mathrm{T}} = \\mathbf{I}$.\n",
        "\n",
        "One of the best-known algebraic feature construction methods is principal component analysis (PCA) - an unsupervised method for dimensionality reduction. Principal _components_ are new features constructed as linear combinations of the given features. The first principal component is given by the direction of maximum variance in the data; the second principal component is the direction of maximum variance orthogonal to the first component, and so on. The _scores_ are numbers (\"scalars\") with which these components can be multiplied to reconstruct the original data. By taking only $k$ $d$-dimensional component vectors and an $n\\times k$ score matrix we have a compressed _approximate_ representation of an $n\\times d$-sized dataset. The information density of all involved numbers/vectors is hopefully higher."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K-rb8nH3wSny"
      },
      "source": [
        "## Load simple numerical dataset\n",
        "\n",
        "Our simple, first dataset consists of four instances of measurements of two features. In terms of the mathematical symbols: this is our $n \\times d = 4 \\times 2$-sized matrix $\\mathbf{X}$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ucRE1KSpwSny"
      },
      "outputs": [],
      "source": [
        "X = np.array([[0.0, 0.4], \n",
        "              [1.0, 2.0], \n",
        "              [2.0, 3.2], \n",
        "              [3.0, 5.3]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RranDBzrwSnz"
      },
      "source": [
        "\n",
        "## Center and scale\n",
        "\n",
        "We need to normalize our dataset if we want PCA to consider all features potentially equally important. Specifically, we want all $d$ features to be centered on zero and have the same standard deviation."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rmli24BGg5iR"
      },
      "source": [
        "### Exercise 1: implement `(un)center_and_(un)scale`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oGHsooLvwSn0"
      },
      "outputs": [],
      "source": [
        "def center_and_scale(X):\n",
        "    \"\"\"`X` is a numpy.ndarray. Rows represent instances. Columns represent feature values.\n",
        "    First output is a modified copy of `X` where mean of individual features is 0 and \n",
        "    standard deviation is 1. The second and third output are original values of said mean and \n",
        "    standard deviation respectively.\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "\n",
        "def uncenter_and_unscale(X_norm, orig_center, orig_std):\n",
        "    \"\"\"`X_norm` is a numpy.ndarray. Rows represent instances. Columns represent feature values,\n",
        "    now with mean 0 and standard deviation 1. Output is a modified copy of `X_norm` with values of\n",
        "    `orig_center` as column means and values of `orig_std` as column standard deviations. (This\n",
        "    function thus reverses `center_and`scale`.)\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "\n",
        "X_norm, orig_center, orig_std = center_and_scale(X)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uv1YeSa4n53r"
      },
      "outputs": [],
      "source": [
        "def test_center_and_scale():\n",
        "    # And uncenter_and_unscale...\n",
        "    X = np.array([[0.0, 0.5, 7], [1.0, 1.5, -10], [2.0, 4.5, 0], [3.0, 5.5, 999]])\n",
        "    X_norm, orig_center, orig_std = center_and_scale(X)\n",
        "    assert helpful_eq(X_norm.mean(), 0)\n",
        "    assert helpful_eq(X_norm.var(), 1)\n",
        "    sk_X_norm = sk_scale(X)\n",
        "    assert helpful_eq(X_norm, sk_X_norm)\n",
        "    assert helpful_eq(orig_center, [1.5, 3.0, 249.])\n",
        "    assert helpful_eq(orig_std, [1.118, 2.062, 433.055])\n",
        "    X = np.dot(np.random.rand(2, 2), np.random.randn(2, 100)).T\n",
        "    csX = center_and_scale(X)\n",
        "    assert helpful_eq(csX[0], sk_scale(X))\n",
        "    assert helpful_eq(uncenter_and_unscale(*csX), X)\n",
        "\n",
        "test_center_and_scale()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S8b1PgAewSn5"
      },
      "source": [
        "## Visualize our goal\n",
        "\n",
        "We are using https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html. You can use this cell as inspiration for the other programming exercises if you didn't manage to finish the earlier ones."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zHGGtCELwSn5",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "k = 1\n",
        "sk_pca = PCA(n_components=k)\n",
        "sk_X_norm = sk_scale(X)\n",
        "sk_score = sk_pca.fit_transform(sk_X_norm)\n",
        "sk_components = sk_pca.components_\n",
        "X_norm_reconstr = sk_score @ sk_components  # @ is a standard matrix multiplication\n",
        "\n",
        "print(\"Component(s) matrix:\")\n",
        "print(sk_components)\n",
        "print(\"Score matrix\")\n",
        "print(sk_score)\n",
        "print()\n",
        "print(sk_X_norm.round(3))\n",
        "print(f\"Recreation of [centered and scaled matrix above] using {k} component(s):\")\n",
        "print(X_norm_reconstr.round(3))\n",
        "mean_squared_error(sk_X_norm, X_norm_reconstr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oFOJUNrmymPX"
      },
      "source": [
        "\n",
        "### Question 1: why can the two-dimensional data we stored in `X` be reconstructed so well using only one principal component?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "usWY5XLAn53y"
      },
      "source": [
        "## Singular value decomposition\n",
        "\n",
        "PCA can be explained in a number of different ways: here, we will derive it by means of the singular value decomposition (SVD). Any $n$-by-$d$ matrix can be uniquely written as a product of three matrices with special properties: \n",
        "\n",
        "\\begin{equation}\n",
        "\\mathbf{X} = \\mathbf{U} \\mathbf{\\Sigma} \\mathbf{V}^{\\mathrm{T}}\n",
        "\\end{equation}\n",
        "\n",
        "Here, $\\mathbf{U}$ is an $n$-by-$k$ matrix, $\\mathbf{\\Sigma}$ is an $k$-by-$k$ matrix and $\\mathbf{V}$ is an $d$-by-$k$ matrix. We will assume that $k = d < n$, which is practically fine but really not in line with the classical definition of an SVD. The NumPy function (https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html) works with a similar assumption/simplification, but goes further by returning only the diagonal elements of $\\mathbf{\\Sigma}$. \n",
        "\n",
        "Furthermore, $\\mathbf{U}$ and $\\mathbf{V}$ are orthogonal (hence rotations) and $\\mathbf{\\Sigma}$ is diagonal (hence a scaling). The columns of $\\mathbf{U}$ and $\\mathbf{V}$ are known as the left and right singular vectors, respectively; and the values on the diagonal of $\\mathbf{\\Sigma}$ are the corresponding singular values.\n",
        "\n",
        "The `R_squared` variable contains the amount of variance the corresponding principal component \"explains\".\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5516tMyZn53z"
      },
      "outputs": [],
      "source": [
        "U, Sigma_elements, V_transpose = np.linalg.svd(X_norm, full_matrices=False)\n",
        "Sigma = np.zeros((U.shape[1], V_transpose.shape[0]))\n",
        "Sigma[:Sigma_elements.size, :Sigma_elements.size] = np.diag(Sigma_elements)\n",
        "S_squared_sum = sum(S_i ** 2 for S_i in Sigma_elements)\n",
        "R_squared = [S_i ** 2 / S_squared_sum for S_i in Sigma_elements]\n",
        "print(U)\n",
        "print(Sigma)\n",
        "print(Sigma_elements)\n",
        "print(V_transpose)\n",
        "print(R_squared)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZHPuecBDwSn3"
      },
      "source": [
        "## Obtain $k$ principal components and score matrix\n",
        "\n",
        "We will name the $k\\times d$ principal components matrix \"components\" and the $n\\times k$ score matrix \"score\" so that we don't keep working with unidiomatic variables. The components matrix can be multiplied with the score matrix to create a reconstruction of $\\mathbf{X}$. Here's an even bigger hint: the score matrix can be found in $\\mathbf{U}\\mathbf{\\Sigma}$ and the components are somewhere in $\\mathbf{V}^{\\mathrm{T}}$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x4i4nhTMhel5"
      },
      "source": [
        "### Exercise 2: implement `svd_to_components_and_score`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xkJV0UqAwSn3"
      },
      "outputs": [],
      "source": [
        "def svd_to_components_and_score(U, Sigma_elements, V_transpose, k):\n",
        "    \"\"\" Input is output of `numpy.linalg.svd(..., full_matrices=False)` and `k` indicating # of \n",
        "    components. First output is `k`-by-d component numpy.ndarray; second output is n-by-`k` score \n",
        "    numpy.ndarray.\n",
        "    \"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "components, score = svd_to_components_and_score(U, Sigma_elements, V_transpose, k)\n",
        "print(\"Component(s) matrix:\")\n",
        "print(components)\n",
        "print(\"Score matrix\")\n",
        "print(score)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Rb7-isB7ymPa"
      },
      "outputs": [],
      "source": [
        "def test_svd_to_components_and_score():\n",
        "    for k in [1, 2]:\n",
        "        X = sk_scale(np.array([[0.0, 0.4],\n",
        "                               [1.0, 1.6],\n",
        "                               [2.0, 2.4],\n",
        "                               [3.0, 3.6]]))\n",
        "        U, S, V_transpose = np.linalg.svd(X, full_matrices=False)\n",
        "        components, score = svd_to_components_and_score(U, S, V_transpose, k)\n",
        "        assert helpful_eq(components[0, 0], 0, failing_is_good=True)\n",
        "        assert helpful_eq(components[0, 0], components[0, 1])\n",
        "        assert helpful_eq(score[0, 0], -score[3, 0])\n",
        "        assert helpful_eq(score[1, 0], -score[2, 0])\n",
        "        np.random.seed(39)\n",
        "        X = sk_scale(np.dot(np.random.rand(2, 2), np.random.randn(2, 100)).T)\n",
        "        sk_pca = PCA(n_components=k)\n",
        "        sk_score = sk_pca.fit_transform(X)\n",
        "        sk_components = sk_pca.components_\n",
        "        U, S, V_transpose = np.linalg.svd(X, full_matrices=False)\n",
        "        components, score = svd_to_components_and_score(U, S, V_transpose, k)\n",
        "        assert helpful_eq(sk_components, components)\n",
        "        assert helpful_eq(sk_score, score)\n",
        "\n",
        "\n",
        "test_svd_to_components_and_score()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XwhmjnYXwSn4"
      },
      "source": [
        "## Reconstruct our data using $k=1$ principal components\n",
        "\n",
        "If you don't have `orig_center` and `orig_std`, you can compare `X_norm` with its reconstruction."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iVgVzucewSn4",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "print(X)\n",
        "print(f\"Recreation of [matrix above] using {k} components:\")\n",
        "X_norm_reconstr = score @ components\n",
        "X_reconstr = uncenter_and_unscale(X_norm_reconstr, orig_center, orig_std)\n",
        "print(X_reconstr.round(3))\n",
        "mean_squared_error(X_norm, X_norm_reconstr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hP9TzzY7Fsr4"
      },
      "source": [
        "### Question 3: make a scatter plot that contains true and reconstructed values $(k=1)$, as well as the single principal component.\n",
        "\n",
        "It is easiest to use the centered and scaled versions of \n",
        "$\\mathbf{X}$ and its reconstruction. You can optionally scale the component visualization with `R_squared` as follows. That is not so informative with $k=1$ but it might be nice to try with $k=2$. Here is some code that you can probably use to get an idea of how:\n",
        "\n",
        "```\n",
        "for i in range(k):\n",
        "    horizontal = R_squared[i] * components[i, 0]\n",
        "    vertical = R_squared[i] * components[i, 1]\n",
        "    plt.arrow(0, 0, horizontal, vertical, head_width=0.05, color=\"red\")\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Av-I5E4lFuNn"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4VRmTKPon54C"
      },
      "source": [
        "## Score and reconstruct _new_ data\n",
        "\n",
        "We will now represent new data in terms of $k$ principal components. For that, we first need to normalize the \"new\" data's feature values using the \"old\" `orig_center` and `orig_std`, so that everything is comparable.\n",
        "\n",
        "We get the new $\\mathbf{V}^{\\mathrm{T}}$ with $\\mathbf{X} (\\mathbf{U}\\mathbf{\\Sigma})^{-1}=\\mathbf{X} (\\mathbf{U}\\mathbf{\\Sigma})^{\\mathrm{T}}$. It could be useful to try to figure out why."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BqDhE13Sh0Dy"
      },
      "source": [
        "### Exercise 3: implement `score_new_data`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VqszvTPQn54D"
      },
      "outputs": [],
      "source": [
        "def score_new_data(new_X, orig_center, orig_std, components):\n",
        "    \"\"\"`new_X` is a numpy.ndarray. Rows represent instances. Columns represent feature values.\n",
        "    First, create a modified copy of `new_X` as if normalizing it using `orig_center` and \n",
        "    `orig_std`. These params refer to feature means and standard deviations respectively. \n",
        "    Output is a representation of this normalized `new_X` in terms of `components`, i.e., a score \n",
        "    matrix.\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "print(X.round(3))\n",
        "new_X = np.array([[0.5, 0.99],\n",
        "                  [2.5, 4.42]])\n",
        "# If you don't have `orig_center` and `orig_std`, you can come up with some fake values.\n",
        "# (You can find the dimensions of these variables in their test function.)\n",
        "new_score = score_new_data(new_X, orig_center, orig_std, components)\n",
        "new_X_norm_reconstr = new_score @ components\n",
        "print(uncenter_and_unscale(new_X_norm_reconstr, orig_center, orig_std).round(3))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "N-WsDaYeymPg"
      },
      "outputs": [],
      "source": [
        "def test_score_new_data():\n",
        "    for k in [1,2,3]:\n",
        "\n",
        "        X = sk_scale(np.array([[0.0, 0.4, -1, 1],\n",
        "                               [1.0, 1.6, -2, 2],\n",
        "                               [2.0, 3.4, -3, 3.5],\n",
        "                               [3.0, 4.6, -4, 5]]))\n",
        "        new_X = sk_scale(np.array([[0.5, 0.99, 1, -2],\n",
        "                                   [2.5, 3.99, -42, -9]]))\n",
        "        sk_pca = PCA(n_components=k)\n",
        "        score = sk_pca.fit_transform(X)\n",
        "        components = sk_pca.components_\n",
        "        new_score = score_new_data(new_X, [1.42, 2.42, 4, 4], [1.11, 1.66, 2, 3], components)\n",
        "        new_X_reconstr = new_score @ components\n",
        "        if k == 1:\n",
        "            assert helpful_eq(new_X_reconstr, \n",
        "                              np.array([[-0.93533774, -0.93443902,  0.93533774, -0.93414063],\n",
        "                                        [-0.09948932, -0.09939373,  0.09948932, -0.09936199]]))\n",
        "        elif k == 2:\n",
        "            assert helpful_eq(new_X_reconstr, \n",
        "                              np.array([[-0.57078589, -1.07598312,  0.57078589, -1.52258937],\n",
        "                                        [ 0.98203221, -0.51931478, -0.98203221, -1.84512183]]))\n",
        "        elif k == 3:\n",
        "            assert helpful_eq(new_X_reconstr, \n",
        "                              np.array([[-0.34009009, -2.06024096,  0.34009009, -1.        ],\n",
        "                                        [ 1.06081081, -0.85542169, -1.06081081, -1.66666667]]))\n",
        "\n",
        "\n",
        "test_score_new_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S_9J_1BOwSn6"
      },
      "source": [
        "## Wrap our work in a `pca` function\n",
        "\n",
        "This is fairly self-explanatory. You could create an object to keep track of so many variables, but you have seen that trick before; a dictionary works too. If you do not break anything that is already here, you could store more relatively small variables in the output."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hAquONvZwSn6"
      },
      "outputs": [],
      "source": [
        "def pca(X, k):\n",
        "    \"\"\"All components are explained above. If you did not complete `center_and_scale` or \n",
        "    `svd_to_components_and_score`, you can fit in the scikit-learn alternatives exemplified \n",
        "    (in the tests) above.\"\"\"\n",
        "    X_norm, orig_center, orig_std = center_and_scale(X)\n",
        "    U, Sigma_elements, V_transpose = np.linalg.svd(X_norm, full_matrices=False)\n",
        "    S_squared_sum = sum(S_i ** 2 for S_i in Sigma_elements)\n",
        "    R_squared = [S_i ** 2 / S_squared_sum for S_i in Sigma_elements]\n",
        "    components, score = svd_to_components_and_score(U, Sigma_elements, V_transpose, k)\n",
        "    return {\n",
        "        'orig_center': orig_center,\n",
        "        'orig_std': orig_std,\n",
        "        'components': components,\n",
        "        'score': score,\n",
        "        'R_squared': np.array(R_squared)\n",
        "    }\n",
        "\n",
        "\n",
        "pca_dict = pca(X, k)\n",
        "for key in pca_dict.keys():\n",
        "    print(key)\n",
        "    print(pca_dict[key].round(4))\n",
        "    print()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "du6bqtJbwSn7"
      },
      "source": [
        "# Big(ger) data\n",
        "\n",
        "We will now use PCA in the pre-processing of a more realistic dataset, and test whether it helps in improving prediction accuracy and/or reducing model fitting and other computational time. We will use the support vector machine classifier that you have seen before. The specific model is not important; that PCA or other pre-processing can help, is. (And, to be fair, PCA and other pre-processing techniques can also have more costs than benefits when using other models.)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nia3F-UCIAxX"
      },
      "source": [
        "## Load/create more complicated dataset\n",
        "\n",
        "There is a little more information about the dataset `load_breast_cancer` loads on https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic) and https://scikit-learn.org/stable/datasets/toy_dataset.html#breast-cancer-dataset. We modify the dataset to consist of more samples through repeats, but add a decent amount of noise (independently sampled from a normal distribution with mean zero and a quarter of that feature's own standard deviation). We also add two categorical features in the columns of `X` that are in `non_num_idx`. \n",
        "\n",
        "So the dataset does not represent what it's based on very closely anymore. But to give a little more information, in `X, y = load_breast_cancer(return_X_y=True)` set,  `X` has 569 instances of 30 numeric features but 10 categories of attributes:\n",
        "\n",
        "0. radius (mean of distances from center to points on the perimeter)\n",
        "1. texture (standard deviation of gray-scale values)\n",
        "2. perimeter\n",
        "3. area\n",
        "4. smoothness (local variation in radius lengths)\n",
        "5. compactness (perimeter^2 / area - 1.0)\n",
        "6. concavity (severity of concave portions of the contour)\n",
        "7. concave points (number of concave portions of the contour)\n",
        "8. symmetry\n",
        "9. fractal dimension (“coastline approximation” - 1)\n",
        "\n",
        "Features are computed from a digitized image of a fine needle aspirate (FNA) (sample) from a tumor. They describe characteristics of the cell nuclei present in the image. The mean, standard error, and mean of the three largest values of these characteristics were computed for each of the 569 images/instances, resulting in 30 features. For instance, field 0 is mean radius, field 10 is radius standard error, field 20 is the mean of the three largest radii.\n",
        "\n",
        "The `y` contains a 0 for malignant samples, and a 1 for benign ones.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cCuVeoBDn54K"
      },
      "outputs": [],
      "source": [
        "def create_big_X_y():\n",
        "    \"\"\"It's not necessary to understand completely what is happening, but notice that \n",
        "    the columns whose indices are in `non_num_idx` represent categorical features.\"\"\"\n",
        "    X, y = load_breast_cancer(return_X_y=True)\n",
        "\n",
        "    # Add some instances\n",
        "    X = np.repeat(X, 10, axis=0)\n",
        "    y = np.repeat(y, 10)\n",
        "\n",
        "    X, y = sk_shuffle(X, y, random_state=11)\n",
        "\n",
        "    # Add two categorical features (that have a 60% prob of being right)\n",
        "    np.random.seed(42)\n",
        "    X = np.hstack((np.random.binomial(1, 0.4 + 0.2 * y).reshape(-1, 1),\n",
        "                   np.random.binomial(1, 0.4 + 0.2 * y).reshape(-1, 1),\n",
        "                   X))\n",
        "    non_num_idx = [0, 1]\n",
        "    \n",
        "    # Add noise to numerical features\n",
        "    num = [i for i in range(X.shape[1]) if i not in non_num_idx]\n",
        "    X[:, num] += np.random.normal(0, 0.25 *\n",
        "                                  X[:, num].std(axis=0), X[:, num].shape)\n",
        "\n",
        "    # Replace some % of numerical values with np.nan\n",
        "    X_num_flat = X[:, num].flatten()\n",
        "    n_to_remove = int(len(X_num_flat) * .1)\n",
        "    to_remove = np.random.permutation(range(len(X_num_flat)))[:n_to_remove]\n",
        "    X_num_flat[to_remove] = np.nan\n",
        "    X[:, num] = X_num_flat.reshape(X[:, num].shape)\n",
        "\n",
        "    return X, y, non_num_idx\n",
        "\n",
        "\n",
        "# big_... name to not confuse with our previous simple `X`.\n",
        "big_X, big_y, non_num_idx = create_big_X_y()\n",
        "big_X.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G2g5edAcMQO8"
      },
      "source": [
        "## Already split off a test set\n",
        "\n",
        "We will perform repeated cross-validation later on, but while we create and test our functions it is already useful to have a separation of training and test set in mind."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v7vSOUxQMRNe",
        "scrolled": false
      },
      "outputs": [],
      "source": [
        "X_train, X_test, y_train, y_test = NFolds(big_X, big_y, seed=5).get_fold(0)\n",
        "print(X_train)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iKgW9ZjrskoK"
      },
      "source": [
        "## Explore the training set"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bSTvz5nOskoK"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gS25qzlWnwrk",
        "pycharm": {
          "name": "#%% md\n"
        }
      },
      "source": [
        "### Question 4: can you come up with a feature transformation (for one or more column(s) of `X`) that would make the principal components express more of the data's variance?\n",
        "\n",
        "Feel free to demonstrate this in a meaningful way."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EL-wAlZ2IEOd"
      },
      "source": [
        "## Replace missing _numerical_ values\n",
        "\n",
        "We will use a simple mean imputation to replace the missing values. Remember that we do not want to use any information from our test set before our final evaluation."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dh0nIGEzyt9R"
      },
      "source": [
        "### Exercise 4: implement `mean_impute`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "D6RCMCqIn54M"
      },
      "outputs": [],
      "source": [
        "def mean_impute(X_train, X_test):\n",
        "    \"\"\"We assume numpy.ndarrays `X_train` and `X_test` only have numpy.nan as missing values in \n",
        "    their numerical features (columns). The output should be modified copies of `X_train` and \n",
        "    `X_test` where the missing values are replaced by the mean of the column of `X_train` in which \n",
        "    the value was missing. (The test might enlighten you.)\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "X_train, X_test = mean_impute(X_train, X_test)\n",
        "print(X_train)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6SurEQy2ymPk"
      },
      "outputs": [],
      "source": [
        "def test_mean_impute():\n",
        "    X_train = np.array([[1, 2, np.nan],\n",
        "                        [3, 4, 5]])\n",
        "    X_test = np.array([[np.nan, 0, 3]])\n",
        "    X_train, X_test = mean_impute(X_train, X_test)\n",
        "    assert helpful_eq(X_train, [[1., 2., 5.], [3., 4., 5.]])\n",
        "    assert helpful_eq(X_test, [[2., 0., 3.]])\n",
        "    big_X, big_y, non_num_idx = create_big_X_y()\n",
        "    X_train, X_test, _, _ = NFolds(big_X, big_y, seed=5).get_fold(0)\n",
        "    X_train, X_test = mean_impute(X_train, X_test)\n",
        "    assert helpful_eq(np.mean(X_train, axis=1)[0], 62.12441680772882)\n",
        "    assert helpful_eq(np.mean(X_test, axis=1)[0], 53.21256512320771)\n",
        "    assert helpful_eq(np.mean(X_train), 57.90520707104668)\n",
        "    assert helpful_eq(np.mean(X_test), 58.20296278479702)\n",
        "\n",
        "test_mean_impute()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X7sn7kqMIJEc"
      },
      "source": [
        "## Choose $k$ components"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b8z_ynjlyy8_"
      },
      "source": [
        "### Exercise 5: implement `only_num`\n",
        "\n",
        "We want to perform PCA on only the numerical part of the dataset."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DM6IiF40aNmT"
      },
      "outputs": [],
      "source": [
        "def only_num(X, non_num_idx):\n",
        "    \"\"\"The columns of numpy.ndarray `X` whose indices are stored in `non_num_idx` are not \n",
        "    numerical. Return a modified copy of `X` that does _not_ contain these columns.\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "print(only_num(X_train, non_num_idx).round(3))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cp5Fgq8VaPw3"
      },
      "outputs": [],
      "source": [
        "def test_only_num():\n",
        "    X = np.array([[-1, 0, 1, 2],\n",
        "                  [-3, 1, 4, 5]])\n",
        "    assert helpful_eq(only_num(X, [1]), np.array([[-1, 1, 2],\n",
        "                                                  [-3, 4, 5]]))\n",
        "    big_X, big_y, non_num_idx = create_big_X_y()\n",
        "    first_num_idx = min([i for i in range(big_X.shape[1]) if i not in non_num_idx])\n",
        "    assert helpful_eq(only_num(big_X, non_num_idx)[0:2, 0], big_X[0:2, first_num_idx])\n",
        "\n",
        "\n",
        "test_only_num()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FBX-zM_nwY8y"
      },
      "source": [
        "### Question 5: what is a good number of principal components to continue with _and why_? (Base your answer only on `X_train`.)\n",
        "\n",
        "A figure may be useful."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dt2uDOTJI0RL"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FU6jqajtQ-Us"
      },
      "outputs": [],
      "source": [
        "k = 2\n",
        "pca_dict = pca(only_num(X_train, non_num_idx), k)\n",
        "for key in pca_dict.keys():\n",
        "    print(key)\n",
        "    print(pca_dict[key].round(4))\n",
        "    print()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9505SsXZIfgX"
      },
      "source": [
        "## Replace `X`'s numerical features with score matrix\n",
        "\n",
        "`X` has many numerical features and a couple of others, whose indices are stored in `non_num_idx`. We only perform PCA on the numerical features. We want to re-group the outcome of the PCA and the non-numerical part of `X` again."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EWSUlbvdzEI4"
      },
      "source": [
        "### Exercise 6: implement `replace_num_with_score`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DgdZe7mDIjAM"
      },
      "outputs": [],
      "source": [
        "def replace_num_with_score(X, non_num_idx, score):\n",
        "    \"\"\"The columns of numpy.ndarray `X` whose indices are stored in `non_num_idx` are not \n",
        "    numerical. Return a modified copy of `X` that has unchanged non-numerical features but in \n",
        "    which the numerical features are replaced by the elements of the `score` numpy.ndarray.\"\"\"\n",
        "    # Assume only first features can be non-numerical to make our life easier\n",
        "    assert all([i == non_num_idx[i] for i in range(len(non_num_idx))])\n",
        "\n",
        "    raise NotImplementedError()\n",
        "\n",
        "print(X_train.round(3))\n",
        "print(pca_dict['score'].round(3))\n",
        "print(replace_num_with_score(X_train, non_num_idx, pca_dict['score']).round(3))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gYBfiJjJymPn"
      },
      "outputs": [],
      "source": [
        "def test_replace_num_with_score():\n",
        "    X = np.array([[1, 2.1, 4.3],\n",
        "                  [0, 4.2, 5.3]])\n",
        "    X = replace_num_with_score(X, [0], [[0,  1], [3,  2]])\n",
        "    assert helpful_eq(X, [[1, 0, 1], [0, 3, 2]])\n",
        "    X = np.array([[1, 2.1, 4.3]])\n",
        "    X = replace_num_with_score(X, [0, 1], [[50]])\n",
        "    assert helpful_eq(X, [[1, 2.1, 50]])\n",
        "\n",
        "\n",
        "test_replace_num_with_score()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HtlLJsden54T"
      },
      "source": [
        "## Wrap pre-processing in a function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "auwFezgnskoP"
      },
      "outputs": [],
      "source": [
        "def preprocess(X_train, X_test, non_num_idx, pca_k=False):\n",
        "    \"\"\"Returns modified copies of `X_train` and `X_test` in which at least the missing values are\n",
        "    imputed. If `pca_k` is not False but a positive integer, perform PCA as above, and replace\n",
        "    numerical features (column indices not in `non_num_idx`) with elements of the score matrix.\n",
        "    The outcome of the `pca` function is returned as third value, or `None` if no PCA is done.\"\"\"\n",
        "    X_train, X_test = mean_impute(X_train, X_test)\n",
        "    if pca_k > 0:\n",
        "        pca_dict = pca(only_num(X_train, non_num_idx), pca_k)\n",
        "        X_train_with_score = replace_num_with_score(X_train, non_num_idx, pca_dict['score'])\n",
        "        X_test_s = score_new_data(only_num(X_test, non_num_idx), pca_dict['orig_center'], \n",
        "                                  pca_dict['orig_std'], pca_dict['components'])\n",
        "        X_test_with_score = replace_num_with_score(X_test, non_num_idx, X_test_s)\n",
        "        return X_train_with_score, X_test_with_score, pca_dict\n",
        "    else:\n",
        "        return X_train, X_test, None"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YInpizUMskoQ"
      },
      "source": [
        "## Evaluate PCA as part of pre-processing\n",
        "\n",
        "We will use repeated cross validation with https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html to fit and predict exactly the same data. By \"repeated\", we mean dividing the data `n_repeats`times into `n_folds` and having each of these folds as test set once. Compared to doing cross validation just once, this gives a stabler estimate of computation time and accuracy/loss. We are interested in estimating and comparing both of these quantities. To do that fairly, we use the same data split and model fitting seeds.\n",
        "\n",
        "We will use the IPython magic command `%%time` in the beginning of a cell to see the time it -- the cell -- takes to execute."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XYBP-XH8skoQ"
      },
      "outputs": [],
      "source": [
        "n_repeats = n_folds = 4  # you may want to lower this to, e.g., 2 while working on the function"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nxHpSMj8zpmu"
      },
      "source": [
        "### Exercise 7: implement `repeated_cross_validation`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JZIzslpOO-bx"
      },
      "outputs": [],
      "source": [
        "def repeated_cross_validation(X, y, n_repeats, n_folds, non_num_idx, pca_k=False):\n",
        "    \"\"\"For i = 0, ...`n_repeats`, create `X` and `y` folds with seed i. Then perform regular\n",
        "    cross validation per set of folds. This regular CV includes \n",
        "    `preprocess(X_train, X_test, non_num_idx, pca_k)` and training an\n",
        "    `SVC(random_state=(i+1)*(j+1))` where `j` is the index of the test fold. Return an \n",
        "    `n_repeats`-by-`n_folds` numpy.ndarray of mean accuracies per test fold.\"\"\"\n",
        "    raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DmbrlgJ9ymPm"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "k = 5\n",
        "accuracies_pca = repeated_cross_validation(big_X, big_y, n_repeats, n_folds, non_num_idx, pca_k=k)\n",
        "print(accuracies_pca.round(3))\n",
        "print(f'Accuracy with PCA: {round(np.mean(accuracies_pca), 3)}' +\n",
        "      f'+- {round(np.std(accuracies_pca), 3)}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "s5joVCisIdRf"
      },
      "outputs": [],
      "source": [
        "def test_repeated_cross_validation():\n",
        "    X, y = load_breast_cancer(return_X_y=True)\n",
        "    test_mean_impute()\n",
        "    nn_idx = []\n",
        "    accuracies = repeated_cross_validation(X, y, 2, 2, nn_idx, pca_k=False)\n",
        "    assert helpful_eq(accuracies, [[0.91549296, 0.89824561],\n",
        "                                   [0.9084507,  0.91578947]])\n",
        "    accuracies_pca = repeated_cross_validation(X, y, 2, 2, nn_idx, pca_k=2)\n",
        "    assert helpful_eq(accuracies_pca, [[0.91549296, 0.92631579],\n",
        "                                       [0.92605634, 0.93684211]])\n",
        "    X, y = load_breast_cancer(return_X_y=True)\n",
        "    test_mean_impute()\n",
        "    nn_idx = [0, 1, 2]\n",
        "    accuracies = repeated_cross_validation(X, y, 4, 3, nn_idx, pca_k=False)\n",
        "    assert helpful_eq(np.array(accuracies).shape, (4, 3))\n",
        "    assert helpful_eq(accuracies, [[0.92592593, 0.89473684, 0.91052632],\n",
        "                                   [0.8994709,  0.89473684, 0.92105263],\n",
        "                                   [0.91005291, 0.89473684, 0.92631579],\n",
        "                                   [0.92592593, 0.91052632, 0.91052632]])\n",
        "    accuracies_pca = repeated_cross_validation(X, y, 4, 3, nn_idx, pca_k=2)\n",
        "    assert helpful_eq(accuracies_pca, [[0.92063492, 0.88947368, 0.91578947],\n",
        "                                       [0.88888889, 0.9,        0.9       ],\n",
        "                                       [0.8994709,  0.92105263, 0.9       ],\n",
        "                                       [0.93121693, 0.87894737, 0.91052632]])\n",
        "\n",
        "\n",
        "test_repeated_cross_validation()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5mMOJ6BIn54U",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "accuracies = repeated_cross_validation(big_X, big_y, n_repeats, n_folds, non_num_idx)\n",
        "print(accuracies.round(3))\n",
        "print(f'Accuracy without PCA: {round(np.mean(accuracies), 3)}' +\n",
        "      f'+- {round(np.std(accuracies), 3)}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XkqSXNyxymPp"
      },
      "source": [
        "### Question 6: how do the computation time and accuracies differ between the repeated evaluations with and without PCA, and how do you explain this?\n",
        "\n",
        "If you use a different $k$ than you chose in question 5, specify it. We only looked at one specific training set in Q5, so your answer there is not necessarily wrong if you pick something else here.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jhKKG0LeAHeu"
      },
      "source": [
        "### Question 7 (bonus): test your theory.\n",
        "\n",
        "You may include _relevant_ (psuedo-)code and figures in your report. It's fine if you (partially) falsify your theory too, as long as it -- the theory -- made sense in the first place."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "asNmOF_FAH4G"
      },
      "source": [
        "### Question 8: can you describe a situation where (or model for which) PCA would not help?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o0VKGRXxLoEF"
      },
      "source": [
        "## Interpret principal components"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "915cw7gkymPp"
      },
      "source": [
        "### Question 9: make scatterplots of 1: (either train, test, or all) points expressed in terms of the first two principal components; 2: the first two numerical features.\n",
        "\n",
        "Make sure the labels of the points are displayed somehow, so that we can get a rough idea of how well the malignant and benign samples can be distinguished using only the first two principal components.\n",
        "\n",
        "This requires you to run PCA with $k \\geq 2$, but you could use https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html on this dataset as well."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-DBIcW-In54W"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yBXDGoMUO-b1"
      },
      "outputs": [],
      "source": [
        "raise NotImplementedError()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GkaLzfDOymPq"
      },
      "outputs": [],
      "source": [
        "for i in range(pca_dict['components'].shape[1]):\n",
        "    print(\"{:02d}\".format(i), pca_dict['components'][:2, i].round(3))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Xhn5VCXwymPq"
      },
      "source": [
        "### Question 10 (2 points): what do the (first two) principal components represent when you think back to what the numerical features are based on?\n",
        "\n",
        "The extra point is for a high-level interpretation of the entire components (the first two respectively)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L1MOASrRLtCv"
      },
      "source": [
        "## Reconstruct numerical features from score matrix\n",
        "\n",
        "Finally, we want to reconstruct our numerical data on its original scale from its score. PCA with $k < d$ is not supposed to be lossless (i.e., perfect) compression, so our results might not look great, but this is a nice exercise.\n",
        "\n",
        "This may seem more daunting than the other programming exercises but it is just a combination of things we have done before. The docstring contains an outline, and the first case in the `test_` function might help as well."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dzk5tb5yD286"
      },
      "source": [
        "### Exercise 8: implement `approx_X_from_X_with_score`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "aK4tRYtNymPo",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "def approx_X_from_X_with_score(X_with_score, pca_dict, non_num_idx):\n",
        "    \"\"\"Obtain the score matrix from the `non_num_idx` columns of `X_with_score`. Then reconstruct\n",
        "    the feature values they originally represented using the PCA outcome stored in `pca_dict`. \n",
        "    Return these values in the same format as the original matrix, without forgetting the \n",
        "    `non_num_idx` columns.\"\"\"\n",
        "    raise NotImplementedError()\n",
        "\n",
        "\n",
        "_, X_test_with_score, pca_dict = preprocess(X_train, X_test, non_num_idx, pca_k=2)\n",
        "X_test_reconstr = approx_X_from_X_with_score(X_test_with_score, pca_dict, non_num_idx)\n",
        "print(X_test[0].round(3))\n",
        "print(X_test_reconstr[0].round(3))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ouzapTMUFG2n"
      },
      "outputs": [],
      "source": [
        "def test_approx_X_from_X_with_score():\n",
        "    X_train = np.array([[0.0, 0.4], \n",
        "                        [1.0, 1.6], \n",
        "                        [2.0, 3.4], \n",
        "                        [3.0, 4.6]])\n",
        "    X_test = np.array([[4, 5.6]])\n",
        "    non_num_idx = []\n",
        "    _, X_test_w_score, pca_dict = preprocess(X_train, X_test, non_num_idx, pca_k=1)\n",
        "    assert helpful_eq(X_test_w_score, [[2.93797197]])\n",
        "    # ^ tests only if mistake is in previous functions\n",
        "    X_test_reconstr = approx_X_from_X_with_score(X_test_w_score, pca_dict, non_num_idx)\n",
        "    assert helpful_eq(X_test_reconstr, [[3.82267078, 5.85623919]])\n",
        "    assert helpful_eq(mean_squared_error(X_test, X_test_reconstr), 0.04855208630760682)\n",
        "\n",
        "    X, y, non_num_idx = create_big_X_y()\n",
        "    X_train, X_test, y_train, y_test = NFolds(X, y, seed=5).get_fold(0)\n",
        "    X_train, X_test = mean_impute(X_train, X_test)\n",
        "    _, X_test_w_score, pca_dict = preprocess(X_train, X_test, non_num_idx, pca_k=2)\n",
        "    X_test_reconstr = approx_X_from_X_with_score(X_test_w_score, pca_dict, non_num_idx)\n",
        "    assert helpful_eq(X_test_reconstr[0:5, -1], \n",
        "                      [0.07272649, 0.08546382, 0.07856379, 0.07790174, 0.10962264])\n",
        "    assert helpful_eq(mean_squared_error(X_test, X_test_reconstr), 2760.755214107711)\n",
        "\n",
        "\n",
        "test_approx_X_from_X_with_score()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9U84291iuN8C"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true
    },
    "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.9.5"
    },
    "toc": {
      "base_numbering": 1,
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {
        "height": "calc(100% - 180px)",
        "left": "10px",
        "top": "150px",
        "width": "250px"
      },
      "toc_section_display": true,
      "toc_window_display": true
    },
    "varInspector": {
      "cols": {
        "lenName": 16,
        "lenType": 16,
        "lenVar": 40
      },
      "kernels_config": {
        "python": {
          "delete_cmd_postfix": "",
          "delete_cmd_prefix": "del ",
          "library": "var_list.py",
          "varRefreshCmd": "print(var_dic_list())"
        },
        "r": {
          "delete_cmd_postfix": ") ",
          "delete_cmd_prefix": "rm(",
          "library": "var_list.r",
          "varRefreshCmd": "cat(var_dic_list()) "
        }
      },
      "position": {
        "height": "358px",
        "left": "1210px",
        "right": "20px",
        "top": "158px",
        "width": "606px"
      },
      "types_to_exclude": [
        "module",
        "function",
        "builtin_function_or_method",
        "instance",
        "_Feature"
      ],
      "window_display": false
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}