{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine Learning and Statistics for Physicists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `np.einsum` to evaluate the tensor expression $g^{il} \\Gamma^m_{ki} x^k$ which arises in [contravariant derivatives in General Relativity](https://en.wikipedia.org/wiki/Christoffel_symbols#Covariant_derivatives_of_tensors).  Note we are using the GR convention that repeated indices (k,l) are summed over."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "67ebba3716136199857aaeefd0595675",
     "grade": false,
     "grade_id": "cell-b10c5a1cfe3128da",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def tensor_expr(g, Gamma, x, D=4):\n",
    "    \"\"\"Evaluate the tensor expression above.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    g : array\n",
    "        Numpy array of shape (D, D)\n",
    "    Gamma : array\n",
    "        Numpy array of shape (D, D, D)\n",
    "    x : array\n",
    "        Numpy array of shape (D,)\n",
    "    D : int\n",
    "        Dimension of input tensors.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    array\n",
    "        Numpy array of shape (D, D) that evaluates the tensor expression.\n",
    "    \"\"\"\n",
    "    assert g.shape == (D, D)\n",
    "    assert Gamma.shape == (D, D, D)\n",
    "    assert x.shape == (D,)\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "f91f3472db3a535b4aa76494682b0bbd",
     "grade": true,
     "grade_id": "cell-dc1412e0ed9e3c8f",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "g = np.arange(4 ** 2).reshape(4, 4)\n",
    "Gamma = np.arange(4 ** 3).reshape(4, 4, 4)\n",
    "x = np.arange(4)\n",
    "y = tensor_expr(g, Gamma, x)\n",
    "assert np.array_equal(\n",
    "    y,\n",
    "    [[ 1680,  3984,  6288,  8592], [ 1940,  4628,  7316, 10004],\n",
    "     [ 2200,  5272,  8344, 11416], [ 2460,  5916,  9372, 12828]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `np.histogram` to calculate the [probability density](https://en.wikipedia.org/wiki/Probability_density_function) that values in an arbitrary input data array fall within user-specified bins. Hint: `np.histogram` does all the work for you with the correct arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0af6aa7170b9f78496e1dbfd925016db",
     "grade": false,
     "grade_id": "cell-366b67031512bdaa",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def estimate_probability_density(data, bins):\n",
    "    \"\"\"Estimate the probability density of arbitrary data.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array\n",
    "        1D numpy array of random values.\n",
    "    bins : array\n",
    "        1D numpy array of N+1 bin edges to use. Must be increasing.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    array\n",
    "        1D numpy array of N probability densities.\n",
    "    \"\"\"\n",
    "    assert np.all(np.diff(bins) > 0)\n",
    "\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "4929cb30d4bdf875f68e0a5bc610159c",
     "grade": true,
     "grade_id": "cell-3add23f80d497553",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "generator = np.random.RandomState(seed=123)\n",
    "data = generator.uniform(size=100)\n",
    "bins = np.linspace(0., 1., 11)\n",
    "rho = estimate_probability_density(data, bins)\n",
    "assert np.allclose(rho, [ 0.6,  0.8,  0.7,  1.7,  1.1,  1.3,  1.6,  0.9,  0.8,  0.5])\n",
    "data = generator.uniform(size=1000)\n",
    "bins = np.linspace(0., 1., 101)\n",
    "rho = estimate_probability_density(data, bins)\n",
    "dx = bins[1] - bins[0]\n",
    "assert np.allclose(dx * rho.sum(), 1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to calculate the [entropy](https://en.wikipedia.org/wiki/Entropy_estimation) $H(\\rho)$ of a binned probability density, defined as:\n",
    "$$\n",
    "H(\\rho) \\equiv -\\sum_i \\rho_i \\log(\\rho_i) \\Delta w_i \\; ,\n",
    "$$\n",
    "where $\\rho_i$ is the binned density in bin $i$ with width $w_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "cdacb2d2924ff83259b567883a4832d0",
     "grade": false,
     "grade_id": "cell-49d830408cabc403",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def binned_entropy(rho, bins):\n",
    "    \"\"\"Calculate the binned entropy.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    rho : array\n",
    "        1D numpy array of densities, e.g., calculated by the previous function.\n",
    "    bins : array\n",
    "        1D numpy array of N+1 bin edges to use. Must be increasing.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float\n",
    "        Value of the binned entropy.\n",
    "    \"\"\"\n",
    "    assert np.all(np.diff(bins) > 0)\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "2de813a17b08d4d867107f8c4d1b1cee",
     "grade": true,
     "grade_id": "cell-7672bc6e182b3f89",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "generator = np.random.RandomState(seed=123)\n",
    "data1 = generator.uniform(size=10000)\n",
    "data2 = generator.uniform(size=10000) ** 4\n",
    "bins = np.linspace(0., 1., 11)\n",
    "rho1 = estimate_probability_density(data1, bins)\n",
    "rho2 = estimate_probability_density(data2, bins)\n",
    "H1 = binned_entropy(rho1, bins)\n",
    "H2 = binned_entropy(rho2, bins)\n",
    "assert np.allclose(H1, -0.000801544)\n",
    "assert np.allclose(H2, -0.699349908)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal (aka Gaussian) distribution is one of the fundamental probability densities that we will encounter often.\n",
    "\n",
    "Implement the function below using `np.random.multivariate_normal` to generate random samples from an arbitrary multidimensional normal distribution, for a specified random seed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "684ae7a3ae31d055a29b739d4739d371",
     "grade": false,
     "grade_id": "cell-7f58b5b612bdc823",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def generate_normal(mu, C, n, seed=123):\n",
    "    \"\"\"Generate random samples from a normal distribution.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    mu : array\n",
    "        1D array of mean values of length N.\n",
    "    C : array\n",
    "        2D array of covariances of shape (N, N). Must be a positive-definite matrix.\n",
    "    n : int\n",
    "        Number of random samples to generate.\n",
    "    seed : int\n",
    "        Random number seed to use.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    array\n",
    "        2D array of shape (n, N) where each row is a random N-dimensional sample.\n",
    "    \"\"\"\n",
    "    assert len(mu.shape) == 1, 'mu must be 1D.'\n",
    "    assert C.shape == (len(mu), len(mu)), 'C must be N x N.'\n",
    "    assert np.all(np.linalg.eigvals(C) > 0), 'C must be positive definite.'\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "21321d2a8fac0ab9bb90bc268fee5ff8",
     "grade": true,
     "grade_id": "cell-f553238e438b8ee6",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "mu = np.array([-1., 0., +1.])\n",
    "C = np.identity(3)\n",
    "C[0, 1] = C[1, 0] = -0.9\n",
    "Xa = generate_normal(mu, C, n=500, seed=1)\n",
    "Xb = generate_normal(mu, C, n=500, seed=1)\n",
    "Xc = generate_normal(mu, C, n=500, seed=2)\n",
    "assert np.array_equal(Xa, Xb)\n",
    "assert not np.array_equal(Xb, Xc)\n",
    "X = generate_normal(mu, C, n=2000, seed=3)\n",
    "assert np.allclose(np.mean(X, axis=0), mu, rtol=0.001, atol=0.1)\n",
    "assert np.allclose(np.cov(X, rowvar=False), C, rtol=0.001, atol=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to visualize a generated 3D dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_3d():\n",
    "    mu = np.array([-1., 0., +1.])\n",
    "    C = np.identity(3)\n",
    "    C[0, 1] = C[1, 0] = -0.9\n",
    "    X = generate_normal(mu, C, n=2000, seed=3)\n",
    "    df = pd.DataFrame(X, columns=('x0', 'x1', 'x2'))\n",
    "    sns.pairplot(df)\n",
    "    \n",
    "visualize_3d()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read about [correlation and covariance](https://en.wikipedia.org/wiki/Covariance_and_correlation), then implement the function below to create a 2x2 covariance matrix given values of $\\sigma_x$, $\\sigma_y$ and the correlation coefficient $\\rho$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "57152504453db2d2284d37b8453a2524",
     "grade": false,
     "grade_id": "cell-04d6325de7d04380",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def create_2d_covariance(sigma_x, sigma_y, rho):\n",
    "    \"\"\"Create and return the 2x2 covariance matrix specified by the input args.\n",
    "    \"\"\"\n",
    "    assert (sigma_x > 0) and (sigma_y > 0), 'sigmas must be > 0.'\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0df02b718d726b157cb7b4785fce9eb0",
     "grade": true,
     "grade_id": "cell-c55a7fcae406eae4",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "assert np.array_equal(create_2d_covariance(1., 1., 0.0), [[1., 0.], [0., 1.]])\n",
    "assert np.array_equal(create_2d_covariance(2., 1., 0.0), [[4., 0.], [0., 1.]])\n",
    "assert np.array_equal(create_2d_covariance(2., 1., 0.5), [[4., 1.], [1., 1.]])\n",
    "assert np.array_equal(create_2d_covariance(2., 1., -0.5), [[4., -1.], [-1., 1.]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following cell that uses your `create_2d_covariance` and `generate_normal` functions to compare the 2D normal distributions with $\\rho = 0$ (blue), $\\rho = +0.9$ (red) and $\\rho = -0.9$ (green).  (You can ignore the harmless `FutureWarning`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_rhos():\n",
    "    mu = np.zeros(2)\n",
    "    sigma_x, sigma_y = 2., 1.\n",
    "    for rho, cmap in zip((0., +0.9, -0.9), ('Blues', 'Reds', 'Greens')):\n",
    "        C = create_2d_covariance(sigma_x, sigma_y, rho)\n",
    "        X = generate_normal(mu, C, 10000)\n",
    "        sns.kdeplot(X[:, 0], X[:, 1], cmap=cmap)\n",
    "    plt.xlim(-4, +4)\n",
    "    plt.ylim(-2, +2)\n",
    "        \n",
    "compare_rhos()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implement the following neural network layer using pytorch:\n",
    "$$\n",
    "\\mathbf{x}_\\text{out} = \\tanh(W \\mathbf{x}_\\text{in} + \\mathbf{b}) \\; ,\n",
    "$$\n",
    "Note that this layer uses (element-wise) *tanh activation* instead of the *relu activation* from the example in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "c3dc2ca82291410d68a552f59ca401a2",
     "grade": false,
     "grade_id": "cell-d4e91c7aff64912e",
     "locked": false,
     "schema_version": 3,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def xout(W, xin, b):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "f9147341895a6bf0d6a75f95b37abaed",
     "grade": true,
     "grade_id": "cell-9e1efa5f935e18ca",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# A correct solution should pass these tests.\n",
    "W = torch.tensor([[1., 2., 3.], [2., -1., 0.]], requires_grad=True)\n",
    "xin = torch.tensor([0.5, -0.5, 1])\n",
    "b = torch.tensor([-1., 1.])\n",
    "assert torch.allclose(xout(W, xin, b), torch.tensor([0.9051, 0.9866]), rtol=1e-4)"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}