{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation\n", "\n", "### [Neil D. Lawrence](http://inverseprobability.com), University of\n", "\n", "Cambridge\n", "\n", "### 2023-10-17" ], "id": "7d829e7d-db1a-43a1-876e-6fc268b99aff" }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Abstract**: This lecture will introduce the notion of simulation and\n", "review the different types of simulation we might use to represent the\n", "physical world." ], "id": "e63964e8-878b-4692-9858-c9134d4540d9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "$$" ], "id": "e8061a4f-9727-4088-9b42-5b6ca0fc97b1" }, { "cell_type": "markdown", "metadata": {}, "source": [ "::: {.cell .markdown}\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Last lecture Carl Henrik introduced you to some of the challenges of\n", "approximate inference. Including the problem of mathematical\n", "tractability. Before that he introduced you to a particular form of\n", "model, the Gaussian process." ], "id": "37507d44-7aa6-4388-904d-5ca779385939" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "plt.rcParams.update({'font.size': 22})" ], "id": "8840554c-5f25-401b-8865-e9526311af4d" }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ], "id": "ee1ad12e-9f53-4dba-9bd7-da79b42daf6e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## notutils\n", "\n", "\\[edit\\]\n", "\n", "This small package is a helper package for various notebook utilities\n", "used below.\n", "\n", "The software can be installed using" ], "id": "c93a4aa0-1028-4846-b521-27e88ab04321" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install notutils" ], "id": "370de4da-a3e7-4472-bc62-20595d62417e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on GitHub:\n", "\n", "\n", "Once `notutils` is installed, it can be imported in the usual manner." ], "id": "313d105c-4f86-4403-88e2-f7aaa00cfb43" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import notutils" ], "id": "aed4898d-bcd2-46b4-9e26-30aa2d3fe969" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## mlai\n", "\n", "\\[edit\\]\n", "\n", "The `mlai` software is a suite of helper functions for teaching and\n", "demonstrating machine learning algorithms. It was first used in the\n", "Machine Learning and Adaptive Intelligence course in Sheffield in 2013.\n", "\n", "The software can be installed using" ], "id": "c5c582a3-90d5-4950-947f-94198e14dcf3" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install mlai" ], "id": "bd2d03ad-e778-4e76-917b-be2e075bdcb6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "from the command prompt where you can access your python installation.\n", "\n", "The code is also available on GitHub: \n", "\n", "Once `mlai` is installed, it can be imported in the usual manner." ], "id": "95a1c820-a947-43e8-b38a-64a71d9a406c" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai" ], "id": "afcb1b94-746b-45fc-9f37-90458cd8bf8b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Game of Life\n", "\n", "\\[edit\\]\n", "\n", "[John Horton Conway](https://en.wikipedia.org/wiki/John_Horton_Conway)\n", "was a mathematician who developed a game known as the Game of Life. He\n", "died in April 2020, but since he invented the game, he was in effect\n", "‘god’ for this game. But as we will see, just inventing the rules\n", "doesn’t give you omniscience in the game.\n", "\n", "The Game of Life is played on a grid of squares, or pixels. Each pixel\n", "is either on or off. The game has no players, but a set of simple rules\n", "that are followed at each turn the rules are." ], "id": "66cf9b4f-5416-41fc-a6bd-badbd9d29fa6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Life Rules\n", "\n", "\\[edit\\]\n", "\n", "John Conway’s game of life is a cellular automaton where the cells obey\n", "three very simple rules. The cells live on a rectangular grid, so that\n", "each cell has 8 possible neighbors.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "*loneliness*\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "Figure: ‘Death’ through loneliness in Conway’s game of life. If a\n", "cell is surrounded by less than three cells, it ‘dies’ through\n", "loneliness.\n", "\n", "The game proceeds in turns, and at each location in the grid is either\n", "alive or dead. Each turn, a cell counts its neighbors. If there are two\n", "or fewer neighbors, the cell ‘dies’ of ‘loneliness’.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "*overcrowding*\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "Figure: ‘Death’ through overpopulation in Conway’s game of life. If a\n", "cell is surrounded by more than three cells, it ‘dies’ through\n", "loneliness.\n", "\n", "If there are four or more neighbors, the cell ‘dies’ from\n", "‘overcrowding’. If there are three neighbors, the cell persists, or if\n", "it is currently dead, a new cell is born.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "*birth*\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "		\n", "
\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "Figure: Birth in Conway’s life. Any position surrounded by precisely\n", "three live cells will give birth to a new cell at the next turn.\n", "\n", "And that’s it. Those are the simple ‘physical laws’ for Conway’s game.\n", "\n", "The game leads to patterns emerging, some of these patterns are static,\n", "but some oscillate in place, with varying periods. Others oscillate, but\n", "when they complete their cycle they’ve translated to a new location, in\n", "other words they move. In Life the former are known as\n", "[oscillators](https://conwaylife.com/wiki/Oscillator) and the latter as\n", "[spaceships](https://conwaylife.com/wiki/Spaceship)." ], "id": "d09f8b00-a8de-46be-a20c-9596c3df8852" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loafers and Gliders\n", "\n", "\\[edit\\]\n", "\n", "John Horton Conway, as the creator of the game of life, could be seen\n", "somehow as the god of this small universe. He created the rules. The\n", "rules are so simple that in many senses he, and we, are all-knowing in\n", "this space. But despite our knowledge, this world can still ‘surprise’\n", "us. From the simple rules, emergent patterns of behaviour arise. These\n", "include static patterns that don’t change from one turn to the next.\n", "They also include, oscillators, that pulse between different forms\n", "across different periods of time. A particular form of oscillator is\n", "known as a ‘spaceship’, this is one that moves across the board as the\n", "game evolves. One of the simplest and earliest spaceships to be\n", "discovered is known as the glider.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "*Glider (1969)*\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "\n", "\n", "\n", "
\n", "\n", "Figure: *Left* A Glider pattern discovered 1969 by Richard K. Guy.\n", "*Right*. John Horton Conway, creator of *Life* (1937-2020). The glider\n", "is an oscillator that moves diagonally after creation. From the simple\n", "rules of Life it’s not obvious that such an object does exist, until you\n", "do the necessary computation.\n", "\n", "The glider was ‘discovered’ in 1969 by Richard K. Guy. What do we mean\n", "by discovered in this context? Well, as soon as the game of life is\n", "defined, objects such as the glider do somehow exist, but the many\n", "configurations of the game mean that it takes some time for us to see\n", "one and know it exists. This means, that despite being the creator,\n", "Conway, and despite the rules of the game being simple, and despite the\n", "rules being deterministic, we are not ‘omniscient’ in any simplistic\n", "sense. It requires computation to ‘discover’ what can exist in this\n", "universe once it’s been defined.\n", "\n", "\n", "\n", "Figure: The Gosper glider gun is a configuration that creates\n", "gliders. A new glider is released after every 30 turns.\n", "\n", "These patterns had to be discovered, in the same way that a scientist\n", "might discover a disease, or an explorer a new land. For example, the\n", "Gosper glider gun was [discovered by Bill Gosper in\n", "1970](https://conwaylife.com/wiki/Bill_Gosper). It is a pattern that\n", "creates a new glider every 30 turns of the game.\n", "\n", "Despite widespread interest in Life, some of its patterns were only very\n", "recently discovered like the Loafer, discovered in 2013 by Josh Ball.\n", "So, despite the game having existed for over forty years, and the rules\n", "of the game being simple, there are emergent behaviors that are unknown.\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "\n", "*Loafer (2013)*\n", "\n", "
\n", "
\n", "\n", "\n", "\n", "
\n", "		\n", "\n", "\n", "\n", "
\n", "\n", "Figure: *Left* A Loafer pattern discovered by Josh Ball in 2013.\n", "*Right*. John Horton Conway, creator of *Life* (1937-2020).\n", "\n", "Once these patterns are discovered, they are combined (or engineered) to\n", "create new Life patterns that do some remarkable things. For example,\n", "there’s a life pattern that runs a Turing machine, or more remarkably\n", "there’s a Life pattern that runs Life itself.\n", "\n", "\n", "\n", "Figure: The Game of Life running in Life. The video is drawing out\n", "recursively showing pixels that are being formed by filling cells with\n", "moving spaceships. Each individual pixel in this game of life is made up\n", "of $2048 \\times 2048$ pixels called an [OTCA\n", "metapixel](https://www.conwaylife.com/wiki/OTCA_metapixel).\n", "\n", "To find out more about the Game of Life you can watch this video by Alan\n", "Zucconi or read his [associated blog\n", "post](https://www.alanzucconi.com/2020/10/13/conways-game-of-life/)." ], "id": "67c944b1-0b28-411c-9c3e-6bbfe0712a2f" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.lib.display import YouTubeVideo\n", "YouTubeVideo('Kk2MH9O4pXY')" ], "id": "1d1576bb-1708-45c7-bb7c-3a04714533ac" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure: An introduction to the Game of Life by Alan Zucconi.\n", "\n", "Contrast this with our situation where in ‘real life’ we don’t know the\n", "simple rules of the game, the state space is larger, and emergent\n", "behaviors (hurricanes, earthquakes, volcanos, climate change) have\n", "direct consequences for our daily lives, and we understand why the\n", "process of ‘understanding’ the physical world is so difficult. We also\n", "see immediately how much easier we might expect the physical sciences to\n", "be than the social sciences, where the emergent behaviors are contingent\n", "on highly complex human interactions." ], "id": "e4e5cd36-73df-4bff-a5f5-e75f1f2ecb59" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Packing Problems\n", "\n", "\\[edit\\]\n", "\n", "\n", "\n", "Figure: Packing 9 squares into a single square. This example is\n", "trivially solved. Credit \n", "\n", "\n", "\n", "Figure: Packing 17 squares into a single square. The optimal solution\n", "is sometimes hard to find. Here the side length of the smallest square\n", "that holds 17 similarly shaped squares is at least 4.675 times the\n", "smaller square. This solution found by John Bidwell in 1997. Credit\n", "\n", "\n", "Another example of a problem where the “physics” is understood because\n", "it’s really mathematics, is packing problems. Here the mathematics is\n", "just geometry, but still we need some form of compute to solve these\n", "problems. [Erich Friedman’s\n", "website](https://erich-friedman.github.io/packing/) contains a host of\n", "these problems, only some of which are analytically tractable.\n", "\n", "\n", "\n", "Figure: Packing 10 squares into a single square. This example is\n", "proven by Walter Stromquist (Stromquist, 1984). Here\n", "$s=3+\\frac{1}{\\sqrt{2}}$. Credit\n", "" ], "id": "91fe7b07-e944-4cee-b01c-bf2143c3516e" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bayesian Inference by Rejection Sampling\n", "\n", "\\[edit\\]\n", "\n", "One view of Bayesian inference is to assume we are given a mechanism for\n", "generating samples, where we assume that mechanism is representing an\n", "accurate view on the way we believe the world works.\n", "\n", "This mechanism is known as our *prior* belief.\n", "\n", "We combine our prior belief with our observations of the real world by\n", "discarding all those prior samples that are inconsistent with our\n", "observations. The *likelihood* defines mathematically what we mean by\n", "inconsistent with the observations. The higher the noise level in the\n", "likelihood, the looser the notion of consistent.\n", "\n", "The samples that remain are samples from the *posterior*.\n", "\n", "This approach to Bayesian inference is closely related to two sampling\n", "techniques known as *rejection sampling* and *importance sampling*. It\n", "is realized in practice in an approach known as *approximate Bayesian\n", "computation* (ABC) or likelihood-free inference.\n", "\n", "In practice, the algorithm is often too slow to be practical, because\n", "most samples will be inconsistent with the observations and as a result\n", "the mechanism must be operated many times to obtain a few posterior\n", "samples.\n", "\n", "However, in the Gaussian process case, when the likelihood also assumes\n", "Gaussian noise, we can operate this mechanism mathematically, and obtain\n", "the posterior density *analytically*. This is the benefit of Gaussian\n", "processes.\n", "\n", "First, we will load in two python functions for computing the covariance\n", "function." ], "id": "7cc0b41e-4ff3-4866-a50c-03562e86336d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai" ], "id": "38927deb-b0dd-4647-a7ae-819e5d94290a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%load -n mlai.Kernel" ], "id": "91256369-cdbf-4a2a-bab9-463aee04c7ee" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load -n mlai.Kernel\n", "class Kernel():\n", " \"\"\"Covariance function\n", " :param function: covariance function\n", " :type function: function\n", " :param name: name of covariance function\n", " :type name: string\n", " :param shortname: abbreviated name of covariance function\n", " :type shortname: string\n", " :param formula: latex formula of covariance function\n", " :type formula: string\n", " :param function: covariance function\n", " :type function: function\n", " :param \\**kwargs:\n", " See below\n", "\n", " :Keyword Arguments:\n", " * \"\"\"\n", "\n", " def __init__(self, function, name=None, shortname=None, formula=None, **kwargs): \n", " self.function=function\n", " self.formula = formula\n", " self.name = name\n", " self.shortname = shortname\n", " self.parameters=kwargs\n", " \n", " def K(self, X, X2=None):\n", " \"\"\"Compute the full covariance function given a kernel function for two data points.\"\"\"\n", " if X2 is None:\n", " X2 = X\n", " K = np.zeros((X.shape[0], X2.shape[0]))\n", " for i in np.arange(X.shape[0]):\n", " for j in np.arange(X2.shape[0]):\n", " K[i, j] = self.function(X[i, :], X2[j, :], **self.parameters)\n", "\n", " return K\n", "\n", " def diag(self, X):\n", " \"\"\"Compute the diagonal of the covariance function\"\"\"\n", " diagK = np.zeros((X.shape[0], 1))\n", " for i in range(X.shape[0]): \n", " diagK[i] = self.function(X[i, :], X[i, :], **self.parameters)\n", " return diagK\n", "\n", " def _repr_html_(self):\n", " raise NotImplementedError" ], "id": "a34fca2d-73c0-4266-b8b3-d1b346ffc70d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mlai" ], "id": "40a86e57-7a61-48b5-abfd-36a6bbf021a6" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%load -n mlai.eq_cov" ], "id": "d27e5278-6733-4965-aa7b-2e01413148f6" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load -n mlai.eq_cov\n", "def eq_cov(x, x_prime, variance=1., lengthscale=1.):\n", " \"\"\"Exponentiated quadratic covariance function.\"\"\"\n", " diffx = x - x_prime\n", " return variance*np.exp(-0.5*np.dot(diffx, diffx)/lengthscale**2)" ], "id": "3b22db63-6c05-42bd-8db3-4587ed106e0c" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kernel = Kernel(function=eq_cov,\n", " name='Exponentiated Quadratic',\n", " shortname='eq', \n", " lengthscale=0.25)" ], "id": "79339d08-eb69-47a8-9e00-2daeba28a970" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we sample from a multivariate normal density (a multivariate\n", "Gaussian), using the covariance function as the covariance matrix." ], "id": "68c700d9-5002-445d-bd3a-47294fcac6ed" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "np.random.seed(10)\n", "import mlai.plot as plot" ], "id": "1e19b848-2012-4a99-b005-51fd6c6f5a47" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot.rejection_samples(kernel=kernel, \n", " diagrams='./gp')" ], "id": "66a4100c-289d-420c-bdaa-48be807852e9" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import notutils as nu\n", "from ipywidgets import IntSlider" ], "id": "91a15118-51a4-4936-a959-aaeb1f2d6a78" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nu.display_plots('gp_rejection_sample{sample:0>3}.png', \n", " directory='./gp', \n", " sample=IntSlider(1,1,5,1))" ], "id": "e3c01a6c-df3a-4909-895a-63e5b23fea11" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "Figure: One view of Bayesian inference is we have a machine for\n", "generating samples (the *prior*), and we discard all samples\n", "inconsistent with our data, leaving the samples of interest (the\n", "*posterior*). This is a rejection sampling view of Bayesian inference.\n", "The Gaussian process allows us to do this analytically by multiplying\n", "the *prior* by the *likelihood*.\n", "\n", "So, Gaussian processes provide an example of a particular type of model.\n", "Or, scientifically, we can think of such a model as a mathematical\n", "representation of a hypothesis around data. The rejection sampling view\n", "of Bayesian inference can be seen as rejecting portions of that initial\n", "hypothesis that are inconsistent with the data. From a Popperian\n", "perspective, areas of the prior space are falsified by the data, leaving\n", "a posterior space that represents remaining plausible hypotheses.\n", "\n", "The flaw with this point of view is that the initial hypothesis space\n", "was also restricted. It only contained functions where the instantiated\n", "points from the function are jointly Gaussian distributed." ], "id": "52d47ce1-2b42-4b8c-ba86-c85653a1663c" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Universe isn’t as Gaussian as it Was\n", "\n", "\\[edit\\]\n", "\n", "The [Planck space\n", "craft](https://en.wikipedia.org/wiki/Planck_(spacecraft)) was a European\n", "Space Agency space telescope that mapped the cosmic microwave background\n", "(CMB) from 2009 to 2013. The [Cosmic Microwave\n", "Background](https://en.wikipedia.org/wiki/Cosmic_microwave_background)\n", "is the first observable echo we have of the big bang. It dates to\n", "approximately 400,000 years after the big bang, at the time the Universe\n", "was approximately $10^8$ times smaller and the temperature of the\n", "Universe was high, around $3 \\times 10^8$ degrees Kelvin. The Universe\n", "was in the form of a hydrogen plasma. The echo we observe is the moment\n", "when the Universe was cool enough for protons and electrons to combine\n", "to form hydrogen atoms. At this moment, the Universe became transparent\n", "for the first time, and photons could travel through space.\n", "\n", "\n", "\n", "Figure: Artist’s impression of the Planck spacecraft which measured\n", "the Cosmic Microwave Background between 2009 and 2013.\n", "\n", "The objective of the Planck spacecraft was to measure the anisotropy and\n", "statistics of the Cosmic Microwave Background. This was important,\n", "because if the standard model of the Universe is correct the variations\n", "around the very high temperature of the Universe of the CMB should be\n", "distributed according to a Gaussian process.[1] Currently our best\n", "estimates show this to be the case (Elsner et al., 2016, 2015; Jaffe et\n", "al., 1998; Pontzen and Peiris, 2010).\n", "\n", "To the high degree of precision that we could measure with the Planck\n", "space telescope, the CMB appears to be a Gaussian process. The\n", "parameters of its covariance function are given by the fundamental\n", "parameters of the Universe, for example the amount of dark matter and\n", "matter in the Universe.\n", "\n", "\n", "\n", "Figure: The cosmic microwave background is, to a very high degree of\n", "precision, a Gaussian process. The parameters of its covariance function\n", "are given by fundamental parameters of the universe, such as the amount\n", "of dark matter and mass.\n", "\n", "[1] Most of my understanding of this is taken from conversations with\n", "Kyle Cranmer, a physicist who makes extensive use of machine learning\n", "methods in his work. See e.g. Mishra-Sharma and Cranmer (2020) from Kyle\n", "and Siddharth Mishra-Sharma. Of course, any errors in the above text are\n", "mine and do not stem from Kyle." ], "id": "f242d16f-9f05-4c10-9d7e-931ae0a62823" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulating a CMB Map\n", "\n", "The simulation was created by [Boris\n", "Leistedt](https://ixkael.github.io/), see the [original Jupyter notebook\n", "here](https://github.com/ixkael/Prob-tools/blob/master/notebooks/The%20CMB%20as%20a%20Gaussian%20Process.ipynb).\n", "\n", "Here we use that code to simulate our own universe and sample from what\n", "it looks like.\n", "\n", "First, we install some specialist software as well as `matplotlib`,\n", "`scipy`, `numpy` we require\n", "\n", "- `camb`: \n", "- `healpy`: " ], "id": "15135212-7067-4ecb-bcb7-3dea35624039" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install camb" ], "id": "b523a61f-e849-4519-a115-afed068b4e46" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%pip install healpy" ], "id": "9e94ca1c-0709-4147-ab36-dec31fddee48" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%config IPython.matplotlib.backend = 'retina'\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "from matplotlib import rc\n", "from cycler import cycler\n", "\n", "import numpy as np\n", "\n", "rc(\"font\", family=\"serif\", size=14)\n", "rc(\"text\", usetex=False)\n", "matplotlib.rcParams['lines.linewidth'] = 2\n", "matplotlib.rcParams['patch.linewidth'] = 2\n", "matplotlib.rcParams['axes.prop_cycle'] =\\\n", " cycler(\"color\", ['k', 'c', 'm', 'y'])\n", "matplotlib.rcParams['axes.labelsize'] = 16" ], "id": "18208ebe-8a0a-4bc7-99ed-8564f5d0034d" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import healpy as hp\n", "\n", "import camb\n", "from camb import model, initialpower" ], "id": "47d31b9d-8f9a-434c-a229-e68a37b8d96b" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we use the theoretical power spectrum to design the covariance\n", "function." ], "id": "f21b702c-f22d-4267-a967-94a0a1d1498a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nside = 512 # Healpix parameter, giving 12*nside**2 equal-area pixels on the sphere.\n", "lmax = 3*nside # band-limit. Should be 2*nside < lmax < 4*nside to get information content." ], "id": "dd5a29b4-fb13-4025-ab57-75112d004556" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we design our Universe. It is parameterized according to the\n", "[$\\Lambda$CDM model](https://en.wikipedia.org/wiki/Lambda-CDM_model).\n", "The variables are as follows. `H0` is the Hubble parameter (in\n", "Km/s/Mpc). The `ombh2` is Physical Baryon density parameter. The `omch2`\n", "is the physical dark matter density parameter. `mnu` is the sum of the\n", "neutrino masses (in electron Volts). `omk` is the $\\Omega_k$ is the\n", "curvature parameter, which is here set to 0, giving the minimal six\n", "parameter Lambda-CDM model. `tau` is the reionization optical depth.\n", "\n", "Then we set `ns`, the “scalar spectral index”. This was estimated by\n", "Planck to be 0.96. Then there’s `r`, the ratio of the tensor power\n", "spectrum to scalar power spectrum. This has been estimated by Planck to\n", "be under 0.11. Here we set it to zero. These parameters are associated\n", "[with inflation](https://en.wikipedia.org/wiki/Primordial_fluctuations)." ], "id": "b3d34ea4-6c71-465b-aa08-26d9caed345f" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Mostly following http://camb.readthedocs.io/en/latest/CAMBdemo.html with parameters from https://en.wikipedia.org/wiki/Lambda-CDM_model\n", "\n", "pars = camb.CAMBparams()\n", "pars.set_cosmology(H0=67.74, ombh2=0.0223, omch2=0.1188, mnu=0.06, omk=0, tau=0.066)\n", "pars.InitPower.set_params(ns=0.96, r=0)" ], "id": "17b29bcb-73e6-4eb6-aa09-db3e3015fd14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having set the parameters, we now use the python software “Code for\n", "Anisotropies in the Microwave Background” to get the results." ], "id": "1df6fe84-1729-4402-976d-22ec3e2b4cb8" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pars.set_for_lmax(lmax, lens_potential_accuracy=0);\n", "results = camb.get_results(pars)\n", "powers = results.get_cmb_power_spectra(pars)\n", "totCL = powers['total']\n", "unlensedCL = powers['unlensed_scalar']\n", "\n", "ells = np.arange(totCL.shape[0])\n", "Dells = totCL[:, 0]\n", "Cells = Dells * 2*np.pi / ells / (ells + 1) # change of convention to get C_ell\n", "Cells[0:2] = 0" ], "id": "730d8ed3-a521-4d1e-969f-03e184684d61" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cmbmap = hp.synfast(Cells, nside, \n", " lmax=lmax, mmax=None, alm=False, pol=False, \n", " pixwin=False, fwhm=0.0, sigma=None, new=False, verbose=True)" ], "id": "abe3dafe-9fb0-4467-8d4b-d11ce0b747e0" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hp.mollview(cmbmap)\n", "fig = plt.gcf()\n", "mlai.write_figure('mollweide-sample-cmb.png',\n", " directory='./physics/')" ], "id": "573f42bc-8a7a-427e-afd5-50c856aa2959" }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Figure: A simulation of the Cosmic Microwave Background obtained\n", "through sampling from the relevant Gaussian process covariance (in polar\n", "co-ordinates).\n", "\n", "The world we see today, of course, is not a Gaussian process. There are\n", "many discontinuities, for example, in the density of matter, and\n", "therefore in the temperature of the Universe.\n", "\n", "$=f\\Bigg($$\\Bigg)$\n", "\n", "Figure: What we observe today is some non-linear function of the\n", "cosmic microwave background.\n", "\n", "We can think of today’s observed Universe, though, as a being a\n", "consequence of those temperature fluctuations in the CMB. Those\n", "fluctuations are only order $10^{-6}$ of the scale of the overall\n", "temperature of the Universe. But minor fluctuations in that density are\n", "what triggered the pattern of formation of the Galaxies. They determined\n", "how stars formed and created the elements that are the building blocks\n", "of our Earth (Vogelsberger et al., 2020).\n", "\n", "Those cosmological simulations are based on a relatively simple set of\n", "‘rules’ that stem from our understanding of natural laws. These ‘rules’\n", "are mathematical abstractions of the physical world. Representations of\n", "behavior in mathematical form that capture the interaction forces\n", "between particles. The grand aim of physics has been to unify these\n", "rules into a single unifying theory. Popular understanding of this quest\n", "developed because of Stephen Hawking’s book, “[A Brief History of\n", "Time](https://en.wikipedia.org/wiki/A_Brief_History_of_Time)”. The idea\n", "of these laws as ‘ultimate causes’ has given them a pseudo religious\n", "feel, see for example Paul Davies’s book “[The Mind of\n", "God](https://en.wikipedia.org/wiki/The_Mind_of_God)” which comes from a\n", "quotation form Stephen Hawking.\n", "\n", "> If we do discover a theory of everything … it would be the ultimate\n", "> triumph of human reason-for then we would truly know the mind of God\n", ">\n", "> Stephen Hawking in *A Brief History of Time* 1988\n", "\n", "This is an entrancing quote, that seems to work well for selling books\n", "(A Brief History of Time sold over 10 million copies), but as Laplace\n", "has already pointed out to us, the Universe doesn’t work quite so simply\n", "as that. Commonly, God is thought to be omniscient, but having a grand\n", "unifying theory alone doesn’t give us omniscience.\n", "\n", "Laplace’s demon still applies. Even if we had a grand unifying theory,\n", "which encoded “all the forces that set nature in motion” we have an\n", "amount of work left to do in any quest for ‘omniscience’.\n", "\n", "> We may regard the present state of the universe as the effect of its\n", "> past and the cause of its future. An intellect which at a certain\n", "> moment would know all forces that set nature in motion, and all\n", "> positions of all items of which nature is composed, …\n", "\n", "> … if this intellect were also vast enough to submit these data to\n", "> analysis, it would embrace in a single formula the movements of the\n", "> greatest bodies of the universe and those of the tiniest atom; for\n", "> such an intellect nothing would be uncertain and the future just like\n", "> the past would be present before its eyes.\n", ">\n", "> — Pierre Simon Laplace (Laplace, 1814)\n", "\n", "We summarized this notion as $$\n", "\\text{data} + \\text{model} \\stackrel{\\text{compute}}{\\rightarrow} \\text{prediction}\n", "$$ As we pointed out, there is an irony in Laplace’s demon forming the\n", "cornerstone of a movement known as ‘determinism’, because Laplace wrote\n", "about this idea in an essay on probabilities. The more important quote\n", "in the essay was" ], "id": "a8af2434-025e-4537-9c82-f035e97c2f36" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Laplace’s Gremlin\n", "\n", "\\[edit\\]\n", "\n", "> The curve described by a simple molecule of air or vapor is regulated\n", "> in a manner just as certain as the planetary orbits; the only\n", "> difference between them is that which comes from our ignorance.\n", "> Probability is relative, in part to this ignorance, in part to our\n", "> knowledge. We know that of three or greater number of events a single\n", "> one ought to occur; but nothing induces us to believe that one of them\n", "> will occur rather than the others. In this state of indecision it is\n", "> impossible for us to announce their occurrence with certainty. It is,\n", "> however, probable that one of these events, chosen at will, will not\n", "> occur because we see several cases equally possible which exclude its\n", "> occurrence, while only a single one favors it.\n", ">\n", "> — Pierre-Simon Laplace (Laplace, 1814), pg 5\n", "\n", "The representation of ignorance through probability is the true message\n", "of Laplace, I refer to this message as “Laplace’s gremlin”, because it\n", "is the gremlin of uncertainty that interferes with the demon of\n", "determinism to mean that our predictions are not deterministic.\n", "\n", "Our separation of the uncertainty into the data, the model and the\n", "computation give us three domains in which our doubts can creep into our\n", "ability to predict. Over the last three lectures we’ve introduced some\n", "of the basic tools we can use to unpick this uncertainty. You’ve been\n", "introduced to, (or have yow reviewed) *Bayes’ rule*. The rule, which is\n", "a simple consequence of the product rule of probability, is the\n", "foundation of how we update our beliefs in the presence of new\n", "information.\n", "\n", "The real point of Laplace’s essay was that we don’t have access to all\n", "the data, we don’t have access to a complete physical understanding, and\n", "as the example of the Game of Life shows, even if we did have access to\n", "both (as we do for “Conway’s universe”) we still don’t have access to\n", "all the compute that we need to make deterministic predictions. There is\n", "uncertainty in the system which means we can’t make precise predictions.\n", "\n", "Gremlins are imaginary creatures used as an explanation of failure in\n", "aircraft, causing crashes. In that sense the Gremlin represents the\n", "uncertainty that a pilot felt about what might go wrong in a plane which\n", "might be “theoretically sound” but in practice is poorly maintained or\n", "exposed to conditions that take it beyond its design criteria. Laplace’s\n", "gremlin is all the things that your model, data and ability to compute\n", "don’t account for bringing about failures in your ability to predict.\n", "Laplace’s gremlin is the uncertainty in the system.\n", "\n", "\n", "\n", "Figure: Gremlins are seen as the cause of a number of challenges in\n", "this World War II poster.\n", "\n", "Carl Henrik described how a prior probability $p(\\boldsymbol{ \\theta})$\n", "represents our hypothesis about the way the world might behave. This can\n", "be combined with a *likelihood* through the process of multiplication.\n", "Correctly normalized, this gives an updated hypothesis that represents\n", "our *posterior* belief about the model in the light of the data.\n", "\n", "There is a nice symmetry between this approach and how Karl Popper\n", "describes the process of scientific discovery. In *Conjectures and\n", "Refutations* (Popper (1963)), Popper describes the process of scientific\n", "discovery as involving hypothesis and experiment. In our description\n", "hypothesis maps onto the *model*. The model is an abstraction of the\n", "hypothesis, represented for example as a set of mathematical equations,\n", "a computational description, or an analogous system (physical system).\n", "The data is the product of previous experiments, our readings, our\n", "observation of the world around us. We can combine these to make a\n", "prediction about what we might expect the future to hold. Popper’s view\n", "on the philosophy of science was that the prediction should be\n", "falsifiable.\n", "\n", "We can see this process as a spiral driving forward, importantly Popper\n", "relates the relationship between hypothesis (model) and experiment\n", "(predictions) as akin to the relationship between the chicken and the\n", "egg. Which comes first? The answer is that they co-evolve together.\n", "\n", "\n", "\n", "Figure: Experiment, analyze and design is a flywheel of knowledge\n", "that is the dual of the model, data and compute. By running through this\n", "spiral, we refine our hypothesis/model and develop new experiments which\n", "can be analyzed to further refine our hypothesis.\n", "\n", "\n", "\n", "Figure: The sets of different models. There are all the models in the\n", "Universe we might like to work with. Then there are those models that\n", "are computable e.g., by a Turing machine. Then there are those which are\n", "analytical tractable. I.e., where the solution might be found\n", "analytically. Finally, there are Gaussian processes, where the joint\n", "distribution of the states in the model is Gaussian.\n", "\n", "The approach we’ve taken to the model so far has been severely limiting.\n", "By constraining ourselves to models for which the mathematics of\n", "probability is tractable, we severely limit what we can say about the\n", "universe.\n", "\n", "Although Bayes’ rule only implies multiplication of probabilities, to\n", "acquire the posterior we also need to normalize. Very often it is this\n", "normalization step that gets in the way. The normalization step involves\n", "integration over the updated hypothesis space, to ensure the updated\n", "posterior prediction is correct.\n", "\n", "We can map the process of Bayesian inference onto the\n", "$\\text{model} + \\text{data}$ perspective in the following way. We can\n", "see the model as the prior, the data as the likelihood and the\n", "prediction as the posterior[1].\n", "\n", "So, if we think of our model as incorporating what we know about the\n", "physical problem of interest (from Newton, or Bernoulli or Laplace or\n", "Einstein or whoever) and the data as being the observations (e.g., from\n", "Piazzi’s telescope or a particle accelerator) then we can make\n", "predictions about what we might expect to happen in the future by\n", "combining the two. It is *those* predictions that Popper sees as\n", "important in verifying the scientific theory (which is incorporated in\n", "the model).\n", "\n", "But while Gaussian processes are highly flexible non-parametric function\n", "models, they are *not* going to be sufficient to capture the type of\n", "physical processes we might expect to encounter in the real world. To\n", "give a sense, let’s consider a few examples of the phenomena we might\n", "want to capture, either in the scientific world, or in real world\n", "decision making.\n", "\n", "[1] We should be careful about such mappings, this is the one I prefer\n", "to think about because I try to think of my modelling assumptions as\n", "being stored in a probabilistic model, which I see as the prior\n", "distribution over what I expect the data to look like. In many domains\n", "of parametric modelling, however, the prior will be specified over the\n", "parameters of a model. In the Gaussian process formalism we’re using,\n", "this mapping is clearer though. The ‘prior’ is the Gaussian process\n", "prior over functions, the data is the relationship between those\n", "functions and observations we make. This mental model will also suit\n", "what follows in terms of our consideration of simulation. But it would\n", "likely confuse someone who had only come to Bayesian inference through\n", "parametric models such a neural network. Note that even in such models,\n", "there will be a way of writing down the decomposition of the model that\n", "is akin to the above, but it might involve writing down intractable\n", "densities, so it’s often avoided." ], "id": "6e2e7ca3-a989-41cb-9142-33cce46851ab" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Precise Physical Laws\n", "\n", "\\[edit\\]\n", "\n", "We’ve already reviewed the importance of Newton’s laws in forging our\n", "view of science: we mentioned the influence [Christiaan\n", "Huygens’](https://en.wikipedia.org/wiki/Christiaan_Huygens) work on\n", "collisions had on Daniel Bernoulli in forming the kinetic theory of\n", "gases. These ideas inform many of the physical models we have today\n", "around a number of natural phenomena. The MET Office supercomputer in\n", "Exeter spends its mornings computing the weather across the world its\n", "afternoons modelling climate scenarios. It uses the same set of\n", "principles that Newton described, and Bernoulli explored for gases. They\n", "are encoded in the Navier-Stokes equations. Differential equations that\n", "govern the flow of compressible and incompressible fluids. As well as\n", "predicting our weather, these equations are used in fluid dynamics\n", "models to understand the flight of aircraft, the driving characteristics\n", "of racing cars and the efficiency of gas turbine engines.\n", "\n", "This broad class of physical models, or ‘natural laws’ is probably the\n", "closest to what Laplace was referring to in the demon. The search for\n", "unifying physical laws that dictate everything we observe around us has\n", "gone on. Alongside Newton we must mention James Clerk Maxwell, who\n", "unified electricity and magnetism in one set of equations that were\n", "inspired by the work and ideas of Michael Faraday. And still today we\n", "look for unifying equations that bring together in a single mathematical\n", "model the ‘natural laws’ we observe. One equation that for Laplace would\n", "be “all forces that set nature in motion”. We can think of this as our\n", "first time of physical model, a ‘precise model’ of the known laws of our\n", "Universe, a model where we expect that the mapping from the mathematical\n", "abstraction to the physical reality is ‘exact’.[1]\n", "\n", "[1] Unfortunately, I have to use the term ‘exact’ loosely here! For\n", "example, most of these laws treat space/time as a continuum. But in\n", "reality, it is quantised. The smallest length we can define is Planck\n", "length ($1.61 \\times 10^{-35}$), and the the smallest time is Planck\n", "time. So even in this exact world of Maxwell and Newton there is an\n", "abstraction." ], "id": "f61dcebf-c544-4566-9d10-9ca4b00586b0" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Abstraction and Emergent Properties\n", "\n", "\n", "\n", "Figure: A scale of different simulations we might be interested in\n", "when modelling the physical world. The scale is $\\log_{10}$ meters. The\n", "scale reflects something about the level of granularity where we might\n", "choose to know “all positions of all items of which nature is\n", "composed”.\n", "\n", "Unfortunately, even if such an equation were to exist, we would be\n", "unlikely to know “all positions of all items of which nature is\n", "composed”. A good example here is computational systems biology. In that\n", "domain we are interested in understanding the underlying function of the\n", "cell. These systems sit somewhere between the two extremes that Laplace\n", "described: “the movements of the greatest bodies of the universe and\n", "those of the smallest atom”.\n", "\n", "When the smallest atom is considered, we need to introduce uncertainty.\n", "We again turn to a different work of Maxwell, building on Bernoulli’s\n", "kinetic theory of gases we end up with probabilities for representing\n", "the location of the ‘molecules of air’. Instead of a deterministic\n", "location for these particles we represent our belief about their\n", "location in a distribution.\n", "\n", "Computational systems biology is a world of micro-machines, built of\n", "three dimensional foldings of strings of proteins. There are spindles\n", "(stators) and rotors (e.g. [ATP\n", "Synthase](https://en.wikipedia.org/wiki/ATP_synthase)), there are small\n", "copying machines (e.g. [RNA\n", "Polymerase](https://en.wikipedia.org/wiki/RNA_polymerase)) there are\n", "sequence to sequence translators\n", "([Ribosomes](https://en.wikipedia.org/wiki/Ribosome)). The cells store\n", "information in DNA but have an ecosystem of structures and messages\n", "being built and sent in proteins and RNA. Unpicking these structures has\n", "been a major preoccupation of biology. That is knowing where the atoms\n", "of these molecules are in the structure, and how the parts of the\n", "structure move when these small micro-machines are carrying out their\n", "roles.\n", "\n", "We understand most (if not all) of the physical laws that drive the\n", "movements of these molecules, but we don’t understand all the actions of\n", "the cell, nor can we intervene reliably to improve things. So, even in\n", "the case where we have a good understanding of the physical laws,\n", "Laplace’s gremlin emerges in our knowledge of “the positions of all\n", "items of which nature is composed”." ], "id": "2ed788f6-2dad-4141-94bd-96cdbe1e5bd1" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Molecular Dynamics Simulations\n", "\n", "By understanding and simulating the physics, we can recreate operations\n", "that are happening at the level of proteins in the human cell.\n", "[V-ATPase](https://en.wikipedia.org/wiki/V-ATPase) is an enzyme that\n", "pumps protons. But at the microscopic level it’s a small machine. It\n", "produces ATP in response to a proton gradient. A paper in *Science\n", "Advances* (Roh et al., 2020) simulates the functioning of these proteins\n", "that operate across the cell membrane. This makes these proteins\n", "difficult to crystallize, the response to this challenge is to use a\n", "simulation which (somewhat) abstracts the processes. You can also check\n", "this [blog\n", "post](https://www6.slac.stanford.edu/news/2020-10-07-first-detailed-look-how-molecular-ferris-wheel-delivers-protons-cellular-factories)\n", "from the paper’s press release.\n", "\n", "\n", "\n", "Figure: The V-ATPase enzyme pumps proteins across membranes. This\n", "molecular dynamics simulation was published in *Science Advances* (Roh\n", "et al., 2020). The scale is roughly $10^{-8} m$." ], "id": "7760492f-dada-4990-ae61-726660fd5371" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quantum Mechanics\n", "\n", "Alternatively, we can drop down a few scales and consider simulation of\n", "the Schrödinger equation. Intractabilities in the many-electron\n", "Schrödinger equation have been addressed using deep neural networks to\n", "speed up the solution enabling simulation of chemical bonds (Pfau et\n", "al., 2020). The [PR-blog post is also\n", "available](https://deepmind.com/blog/article/FermiNet). The paper uses a\n", "neural network to model the quantum state of a number of electrons.\n", "\n", "\n", "\n", "Figure: The many-electron Schrödinger equation is important in\n", "understanding how Chemical bonds are formed.\n", "\n", "Each of these simulations have the same property of being based on a set\n", "of (physical) rules about how particles interact. But one of the\n", "interesting characteristics of such systems is how the properties of the\n", "system are emergent as the dynamics are allowed to continue.\n", "\n", "These properties cannot be predicted without running the physics, or the\n", "equivalently the equation. Computation is required. And often the amount\n", "of computation that is required is prohibitive." ], "id": "8361f628-cc12-4f09-830c-643bdb114dff" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accelerate Programme\n", "\n", "The Computer Lab is hosting a new initiative, funded by Schmidt Futures,\n", "known as the [Accelerate Programme for Scientific\n", "Discovery](https://acceleratescience.github.io/). The aim is to address\n", "scientific challenges, and accelerate the progress of research, through\n", "using tools in machine learning.\n", "\n", "We now have four fellows appointed, each of whom works at the interface\n", "of machine learning and scientific discovery. They are using the ideas\n", "around machine learning modelling to drive their scientific research.\n", "\n", "For example, [Bingqing\n", "Cheng](https://acceleratescience.github.io/team/bingqing-cheng.html),\n", "one of the Department’s former DECAF Fellows has used neural network\n", "accelerated molecular dynamics simulations to understand a new form of\n", "metallic hydrogen, likely to occur at the heart of stars (Cheng et al.,\n", "2020). The University’s [press release is\n", "here](https://www.cam.ac.uk/research/news/ai-shows-how-hydrogen-becomes-a-metal-inside-giant-planets).\n", "\n", "On her website Bingqing quotes Paul Dirac.\n", "\n", "> The fundamental laws necessary for the mathematical treatment of a\n", "> large part of physics and the whole of chemistry are thus completely\n", "> known, and the difficulty lies only in the fact that application of\n", "> these laws leads to equations that are too complex to be solved.\n", "\n", "> ..approximate practical methods of applying quantum mechanics should\n", "> be developed, which can lead to an explanation of the main features of\n", "> complex atomic systems without too much computation.\n", ">\n", "> — Paul Dirac (6 April 1929)\n", "\n", "Bingqing has now taken a position at IST Austria.\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Vincent Dutordoir\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Sarah Morgan\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Soumya Banerjee\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "Sam Nallaperuma\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Our four current Accelerate fellows are [Challenger\n", "Mishra](https://acceleratescience.github.io/team/challenger-mishra.html),\n", "a physicist interested in string theory and quantizing gravity. [Sarah\n", "Morgan](https://acceleratescience.github.io/team/sarah-morgan.html) from\n", "the Brain Mapping Unit, who is focused on predicting psychosis\n", "trajectories, [Soumya\n", "Bannerjee](https://acceleratescience.github.io/team/soumya-banerjee.html)\n", "who focuses on complex systems and healthcare and [Sam\n", "Nallaperuma](https://acceleratescience.github.io/team/sam-nallaperuma.html)\n", "who the interface of machine learning and biology with particular\n", "interests in emergent behavior in complex systems.\n", "\n", "For those interested in Part III/MPhil projects, you can see their\n", "project suggestions on [this\n", "page](https://mlatcl.github.io/internal/projects/)." ], "id": "463305a7-acc8-4f66-8a98-79eeb345e940" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Related Approaches\n", "\n", "While this module is mainly focusing on emulation as a route to bringing\n", "machine learning closer to the physical world, I don’t want to give the\n", "impression that’s the only approach. It’s worth bearing in mind three\n", "important domains of machine learning (and statistics) that we also\n", "could have explored.\n", "\n", "- Probabilistic Programming\n", "- Approximate Bayesian Computation\n", "- Causal inference\n", "\n", "Each of these domains also brings a lot to the table in terms of\n", "understanding the physical world." ], "id": "406fb23d-dacc-43de-8611-17e2bb8f11cb" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Probabilistic Programming\n", "\n", "Probabilistic programming is an idea that, from our perspective, can be\n", "summarized as follows. What if, when constructing your simulator, or\n", "your model, you used a programming language that was aware of the state\n", "variables and the probability distributions. What if this language could\n", "‘compile’ the program into code that would automatically compute the\n", "Bayesian posterior for you?\n", "\n", "This is the objective of probabilistic programming. The idea is that you\n", "write your model in a language, and that language is automatically\n", "converted into the different modelling codes you need to perform\n", "Bayesian inference.\n", "\n", "The ideas for probabilistic programming originate in\n", "[BUGS](https://www.mrc-bsu.cam.ac.uk/software/bugs/). The software was\n", "developed at the MRC Biostatistics Unit here in Cambridge in the early\n", "1990s, by among others, [Sir David\n", "Spiegelhalter](https://en.wikipedia.org/wiki/David_Spiegelhalter). Carl\n", "Henrik covered in last week’s lecture some of the approaches for\n", "approximate inference. BUGS uses Gibbs sampling. Gibbs sampling,\n", "however, can be slow to converge when there are strong correlations in\n", "the posterior between variables.\n", "\n", "The descendent of BUGS that is probably most similar in the spirit of\n", "its design is [Stan](https://mc-stan.org/). Stan came from researchers\n", "at Columbia University and makes use of a variant of Hamiltonian Monte\n", "Carlo called the No-U-Turn sampler. It builds on automatic\n", "differentiation for the gradients it needs. It’s all written in C++ for\n", "speed, but has interfaces to Python, R, Julia, MATLAB etc. Stan has been\n", "highly successful during the Coronavirus pandemic, with a number of\n", "epidemiological simulations written in the language, for example see\n", "this [blog\n", "post](https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html).\n", "\n", "Other probabilistic programming languages of interest include those that\n", "make use of variational approaches (such as [pyro](https://pyro.ai/))\n", "and allow use of neural network components.\n", "\n", "One important probabilistic programming language being developed is\n", "[Turing](https://turinglang.org/stable/), one of the key developers is\n", "[Hong Ge](https://mlg.eng.cam.ac.uk/hong/) who is a Senior Research\n", "Associate in Cambridge Engineering." ], "id": "4bbd2674-9f01-48af-94b7-ca2383b613ce" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximate Bayesian Computation\n", "\n", "We reintroduced Gaussian processes at the start of this lecture by\n", "sampling from the Gaussian process and matching the samples to data,\n", "discarding those that were distant from our observations. This approach\n", "to Bayesian inference is the starting point for *approximate Bayesian\n", "computation* or ABC.\n", "\n", "The idea is straightforward, if we can measure ‘closeness’ in some\n", "relevant fashion, then we can sample from our simulation, compare our\n", "samples to real world data through ‘closeness measure’ and eliminate\n", "samples that are distant from our data. Through appropriate choice of\n", "closeness measure, our samples can be viewed as coming from an\n", "approximate posterior.\n", "\n", "My Sheffield colleague, Rich Wilkinson, was one of the pioneers of this\n", "approach during his PhD in the Statslab here in Cambridge. You can hear\n", "Rich talking about ABC at NeurIPS in 2013 here." ], "id": "d00f0566-925f-46f1-80e1-518c83cd5d1a" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.lib.display import YouTubeVideo\n", "YouTubeVideo('sssbLkn2JjI')" ], "id": "fb37116f-97e6-45db-b8f7-c94efba58885" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure: Rich Wilkinson giving a Tutorial on ABC at NeurIPS in 2013.\n", "Unfortunately, they’ve not synchronised the slides with the tutorial.\n", "You can find the slides [separately\n", "here](http://media.nips.cc/Conferences/2013/Video/Tutorial2B.pdf)." ], "id": "b62dc665-d31c-4597-9c36-ec7797d0e4f9" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Causality" ], "id": "8b684678-25b2-4164-b7f3-b0083255f793" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.lib.display import YouTubeVideo\n", "YouTubeVideo('yksduYxEusQ')" ], "id": "4d39650b-5572-4447-8e50-355894adf651" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure: Judea Pearl and Elias Bareinboim giving a Tutorial on\n", "Causality at NeurIPS in 2013. Again, the slides aren’t synchronised, but\n", "you can find them separately\n", "[here](http://media.nips.cc/Conferences/2013/nips-dec2013-pearl-bareinboim-tutorial-full.pdf).\n", "\n", "All these approaches offer a lot of promise for developing machine\n", "learning at the interface with science but covering each in detail would\n", "require four separate modules. We’ve chosen to focus on the emulation\n", "approach, for two principal reasons. Firstly, it’s conceptual\n", "simplicity. Our aim is to replace all or part of our simulation with a\n", "machine learning model. Typically, we’re going to want uncertainties as\n", "part of that representation. That explains our focus on Gaussian process\n", "models. Secondly, the emulator method is flexible. Probabilistic\n", "programming requires that the simulator has been built in a particular\n", "way, otherwise we can’t compile the program. Finally, the emulation\n", "approach can be combined with any of the existing simulation approaches.\n", "For example, we might want to write our emulators as probabilistic\n", "programs. Or we might do causal analysis on our emulators, or we could\n", "speed up the simulation in ABC through emulation." ], "id": "dcc3af5a-989e-42da-9ec7-a5c842233619" }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Conclusion\n", "\n", "We’ve introduced the notion of a simulator. A body of computer code that\n", "expresses our understanding of a particular physical system. We\n", "introduced such simulators through *physical laws*, such as laws of\n", "gravitation or electro-magnetism. But we soon saw that in many\n", "simulations those laws become abstracted, and the simulation becomes\n", "more phenomological.\n", "\n", "Even full knowledge of all laws does not give us access to ‘the mind of\n", "God’, because we are lacking information about the data, and we are\n", "missing the compute. These challenges further motivate the need for\n", "abstraction, and we’ve seen examples of where such abstractions are used\n", "in practice.\n", "\n", "The example of Conway’s Game of Life highlights how complex emergent\n", "phenomena can require significant computation to explore." ], "id": "c7ee482a-fca1-4af2-9306-abde36f6424a" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thanks!\n", "\n", "For more information on these subjects and more you might want to check\n", "the following resources.\n", "\n", "- twitter: [@lawrennd](https://twitter.com/lawrennd)\n", "- podcast: [The Talking Machines](http://thetalkingmachines.com)\n", "- newspaper: [Guardian Profile\n", " Page](http://www.theguardian.com/profile/neil-lawrence)\n", "- blog:\n", " [http://inverseprobability.com](http://inverseprobability.com/blog.html)" ], "id": "62423698-4c26-4721-864d-73a31f3a74dd" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ], "id": "13654830-fd31-4959-b105-29e62ce9d395" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cheng, B., Mazzola, G., Pickard, C.J., Ceriotti, M., 2020. Evidence for\n", "supercritical behaviour of high-pressure liquid hydrogen. Nature 585,\n", "217–220. \n", "\n", "Elsner, F., Leistedt, B., Peiris, H.V., 2016. Unbiased pseudo-$C_\\ell$\n", "power spectrum estimation with mode projection. Monthly Notices of the\n", "Royal Astronomical Society 465, 1847–1855.\n", "\n", "\n", "Elsner, F., Leistedt, B., Peiris, H.V., 2015. Unbiased methods for\n", "removing systematics from galaxy clustering measurements. Monthly\n", "Notices of the Royal Astronomical Society 456, 2095–2104.\n", "\n", "\n", "Jaffe, A.H., Bond, J.R., Ferreira, P.G., Knox, L.E., 1998. CMB\n", "likelihood functions for beginners and experts, in: AIP Conf. Proc.\n", "\n", "\n", "Laplace, P.S., 1814. Essai philosophique sur les probabilités, 2nd ed.\n", "Courcier, Paris.\n", "\n", "Mishra-Sharma, S., Cranmer, K., 2020. [Semi-parametric $\\gamma$-ray\n", "modeling with Gaussian processes and variational\n", "inference](https://arxiv.org/abs/2010.10450).\n", "\n", "Pfau, D., Spencer, J.S., Matthews, A.G.D.G., Foulkes, W.M.C., 2020. Ab\n", "initio solution of the many-electron Schrödinger equation with deep\n", "neural networks. Phys. Rev. Research 2, 033429.\n", "\n", "\n", "Pontzen, A., Peiris, H.V., 2010. The cut-sky cosmic microwave background\n", "is not anomalous. Phys. Rev. D 81, 103008.\n", "\n", "\n", "Popper, K.R., 1963. Conjectures and refutations: The growth of\n", "scientific knowledge. Routledge, London.\n", "\n", "Roh, S.-H., Shekhar, M., Pintilie, G., Chipot, C., Wilkens, S.,\n", "Singharoy, A., Chiu, W., 2020. Cryo-EM and MD infer water-mediated\n", "proton transport and autoinhibition mechanisms of Vo complex. Science\n", "Advances 6. \n", "\n", "Stromquist, W.R., 1984. Packing unit squares inside squares, III. Daniel\n", "H. Wagner Associates.\n", "\n", "Vogelsberger, M., Marinacci, F., Torrey, P., Puchwei, E., 2020.\n", "Cosmological simulations of galaxy formation. Nature Reviews Physics 2,\n", "42–66. " ], "id": "5d798d4d-9e76-4fd2-be42-9099a96ed672" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": {} }