diff --git a/LICENSE b/LICENSE index 3604fda..572724f 100644 --- a/LICENSE +++ b/LICENSE @@ -26,4 +26,4 @@ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md index f31fa97..58d7715 100644 --- a/README.md +++ b/README.md @@ -16,17 +16,17 @@

-This is an exemplar project designed to showcase best practices in developing scientific software as part of the ReCoDE Project at Imperial College London. +This is an exemplar project designed to showcase best practices in developing scientific software as part of the ReCoDE Project at Imperial College London. **You do not need to know or care about Markov Chain Monte Carlo for this to be useful to you.** Rather this project is primarily designed to showcase the tools and practices available to you when developing scientific software projects. Maybe you are a PhD student just starting, or a researcher just about to embark on a larger scale software project - there should be something interesting here for you. ## Table of contents -1. [A short introduction][intro] +1. [A short introduction][intro] 1. [Organising code and python packaging][packaging] 1. [Testing your code][testing] -1. Planning the project, MVPs, Premature Optimisation, +1. Planning the project, MVPs, Premature Optimisation, 1. Planning out a larger software project 1. Using Jupyter Notebooks during development 1. Documentation @@ -45,8 +45,8 @@ When you're ready to dive in you have three options: ### 2. Clone the repo and run the jupyter notebooks locally. (Faster but requires you have python/jupyter installed) ``` -git clone -cd +git clone +cd pip install -r requirements.txt jupyter lab ``` @@ -60,7 +60,7 @@ jupyter lab ├── CITATION.cff # This file describes how to cite the work contained in this repository. ├── LICENSE # Outlines what legal rights you have to use this software. ├── README.md # You are here! -├── code +├── code │   ├── README.md # Human readable information about the little python package in here │   ├── pyproject.toml # Machine readable information about that same package │   ├── setup.cfg # Tells Pip how to install this package @@ -71,7 +71,7 @@ jupyter lab ``` ## External Resources -- [The Turing Way](https://the-turing-way.netlify.app/) has tons of great resources on the topics discussed here. +- [The Turing Way](https://the-turing-way.netlify.app/) has tons of great resources on the topics discussed here. - [Intermediate Research Software Development in Python](https://carpentries-incubator.github.io/python-intermediate-development/index.html) diff --git a/code/pyproject.toml b/code/pyproject.toml index fa7093a..b0f0765 100644 --- a/code/pyproject.toml +++ b/code/pyproject.toml @@ -1,3 +1,3 @@ [build-system] requires = ["setuptools>=42"] -build-backend = "setuptools.build_meta" \ No newline at end of file +build-backend = "setuptools.build_meta" diff --git a/code/setup.cfg b/code/setup.cfg index baf5802..845e7b4 100644 --- a/code/setup.cfg +++ b/code/setup.cfg @@ -13,10 +13,10 @@ classifiers = Operating System :: OS Independent [options] -package_dir = +package_dir = = src packages = find: python_requires = >=3.6 [options.packages.find] -where = src \ No newline at end of file +where = src diff --git a/code/src/MCFF/ising_model.py b/code/src/MCFF/ising_model.py index a0ce7e4..e846fd2 100644 --- a/code/src/MCFF/ising_model.py +++ b/code/src/MCFF/ising_model.py @@ -7,15 +7,18 @@ States are represented by NxM numpy arrays with values +/-1. These functions are import numpy as np from numba import jit + def all_up_state(N): "The NxN all up state of the Ising model." return np.ones([N, N]) + def all_down_state(N): "The NxN all down state of the Ising model." return -np.ones([N, N]) -def random_state(N, magnetisation = 0.5, rng = np.random.default_rng()): + +def random_state(N, magnetisation=0.5, rng=np.random.default_rng()): r"""Return the a random NxN state of the Ising model. Parameters @@ -32,13 +35,13 @@ def random_state(N, magnetisation = 0.5, rng = np.random.default_rng()): np.array A random NxN state. """ - return np.random.choice([+1, -1], size = (N,N), p = magnetisation) + return np.random.choice([+1, -1], size=(N, N), p=magnetisation) # nopython=True instructs numba that we want all use of the python interpreter to removed from this function # numba will throw and error if it cannot achieve this. # nogil = True says that numba can release python's Global Interpreter Lock, read more about that here: https://numba.pydata.org/numba-doc/latest/user/jit.html#nogil -@jit(nopython=True, nogil = True) +@jit(nopython=True, nogil=True) def energy(state): r"""Compute the energy of a state of the Ising Model with open boundary conditions. @@ -57,17 +60,18 @@ def energy(state): for j in range(M): # handle the north and south neighbours if 0 <= (i + 1) < N: - E -= state[i,j] * state[i+1, j] - + E -= state[i, j] * state[i + 1, j] + # handle the east and west neighbours if 0 <= (j + 1) < M: - E -= state[i,j] * state[i, j+1] - - return 2*E / (N*M) + E -= state[i, j] * state[i, j + 1] + + return 2 * E / (N * M) + def energy_numpy(state): r"""Compute the energy of a state of the Ising Model with open boundary conditions. - + An alternate implementation of `energy` using Numpy array operations, used for comparison and correctness checks. Parameters @@ -79,6 +83,5 @@ def energy_numpy(state): float64 The interaction energy per site. """ - E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) - return 2*E / np.product(state.shape) - + E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) + return 2 * E / np.product(state.shape) diff --git a/code/src/MCFF/mcmc.py b/code/src/MCFF/mcmc.py index 27ce015..682e57d 100644 --- a/code/src/MCFF/mcmc.py +++ b/code/src/MCFF/mcmc.py @@ -3,7 +3,8 @@ """ -def mcmc(initial_state, steps, T, rng = np.random.default_rng()): + +def mcmc(initial_state, steps, T, rng=np.random.default_rng()): r"""Return the state after performing `steps` monte carlo steps starting from `initial_state` at temperature `T`. Parameters @@ -24,18 +25,18 @@ def mcmc(initial_state, steps, T, rng = np.random.default_rng()): """ N, M = initial_state.shape assert N == M - + current_state = initial_state.copy() E = N**2 * energy(state) for i in range(steps): - i, j = rng.integers(N, size = 2) - + i, j = rng.integers(N, size=2) + new_state = current_state.copy() - new_state[i,j] *= -1 + new_state[i, j] *= -1 new_E = N**2 * energy(new_state) - + if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform(): current_state = new_state E = new_E - - return current_state \ No newline at end of file + + return current_state diff --git a/code/tests/test_energy.py b/code/tests/test_energy.py index 31ba875..9a82cc8 100644 --- a/code/tests/test_energy.py +++ b/code/tests/test_energy.py @@ -6,22 +6,27 @@ from MCFF.ising_model import all_up_state, all_down_state, random_state from MCFF.ising_model import energy, energy_numpy # create some test states -all_up = np.ones([100,100]) -all_down = -np.ones([100,100]) -random = np.random.choice([-1, 1], size = (100,100)) -custom = (1 - 2*imread(Path(__file__).parents[2]/'learning/data/test_state.png')[:, :, 0]) #load a 100x100 png, take the red channel, remap 0,1 to -1,1 +all_up = np.ones([100, 100]) +all_down = -np.ones([100, 100]) +random = np.random.choice([-1, 1], size=(100, 100)) +custom = ( + 1 - 2 * imread(Path(__file__).parents[2] / "learning/data/test_state.png")[:, :, 0] +) # load a 100x100 png, take the red channel, remap 0,1 to -1,1 states = [all_up, all_down, random, custom] + def E_prediction_all_the_same(L): "The exact energy in for the case where all spins are up or down" - return -(4*(L - 2)**2 + 12*(L-2) + 8) / L**2 + return -(4 * (L - 2) ** 2 + 12 * (L - 2) + 8) / L**2 + def test_exact_energies(): for state in [all_up, all_down]: L = state.shape[0] assert energy(state) == E_prediction_all_the_same(L) + def test_energy_implementations(): for state in states: assert np.allclose(energy(state), energy_numpy(state)) diff --git a/code/tests/test_energy_using_hypothesis.py b/code/tests/test_energy_using_hypothesis.py index ecef640..f8e0abc 100644 --- a/code/tests/test_energy_using_hypothesis.py +++ b/code/tests/test_energy_using_hypothesis.py @@ -5,8 +5,13 @@ from hypothesis.extra import numpy as hnp from MCFF.ising_model import energy, energy_numpy -@given(hnp.arrays(dtype = int, - shape = hnp.array_shapes(min_dims = 2, max_dims = 2), - elements = st.sampled_from([1, -1]))) + +@given( + hnp.arrays( + dtype=int, + shape=hnp.array_shapes(min_dims=2, max_dims=2), + elements=st.sampled_from([1, -1]), + ) +) def test_generated_states(state): - assert np.allclose(energy(state), energy_numpy(state)) \ No newline at end of file + assert np.allclose(energy(state), energy_numpy(state)) diff --git a/learning/01 Introduction.ipynb b/learning/01 Introduction.ipynb index 0f6bc98..b73b06d 100644 --- a/learning/01 Introduction.ipynb +++ b/learning/01 Introduction.ipynb @@ -51,9 +51,13 @@ "\n", "# This loads some custom styles for matplotlib\n", "import json, matplotlib\n", - "with open(\"assets/matplotlibrc.json\") as f: matplotlib.rcParams.update(json.load(f))\n", "\n", - "np.random.seed(42) #This makes our random numbers reproducable when the notebook is rerun in order" + "with open(\"assets/matplotlibrc.json\") as f:\n", + " matplotlib.rcParams.update(json.load(f))\n", + "\n", + "np.random.seed(\n", + " 42\n", + ") # This makes our random numbers reproducable when the notebook is rerun in order" ] }, { @@ -84,13 +88,16 @@ } ], "source": [ - "state = np.random.choice([-1, 1], size = (100,100))\n", + "state = np.random.choice([-1, 1], size=(100, 100))\n", + "\n", + "\n", + "def show_state(state, ax=None):\n", + " if ax is None:\n", + " f, ax = plt.subplots()\n", + " ax.matshow(state, cmap=\"Greys\", vmin=-1, vmax=1)\n", + " ax.set(xticks=[], yticks=[])\n", + "\n", "\n", - "def show_state(state, ax = None):\n", - " if ax is None: f, ax = plt.subplots()\n", - " ax.matshow(state, cmap = \"Greys\", vmin = -1, vmax = 1)\n", - " ax.set(xticks = [], yticks = [])\n", - " \n", "show_state(state)" ] }, @@ -132,18 +139,21 @@ } ], "source": [ - "all_up = np.ones([100,100])\n", - "all_down = -np.ones([100,100])\n", - "random = np.random.choice([-1, 1], size = (100,100))\n", + "all_up = np.ones([100, 100])\n", + "all_down = -np.ones([100, 100])\n", + "random = np.random.choice([-1, 1], size=(100, 100))\n", "\n", "from matplotlib.image import imread\n", - "custom = (1 - 2*imread('data/test_state.png')[:, :, 0]) #load a 100x100 png, take the red channel, remap 0,1 to -1,1\n", + "\n", + "custom = (\n", + " 1 - 2 * imread(\"data/test_state.png\")[:, :, 0]\n", + ") # load a 100x100 png, take the red channel, remap 0,1 to -1,1\n", "\n", "states = [all_up, all_down, random, custom]\n", "\n", - "f, axes = plt.subplots(ncols = 4, figsize = (20,5))\n", + "f, axes = plt.subplots(ncols=4, figsize=(20, 5))\n", "for ax, state in zip(axes, states):\n", - " show_state(state, ax = ax)" + " show_state(state, ax=ax)" ] }, { @@ -181,7 +191,8 @@ ], "source": [ "def E_prediction_all_the_same(L):\n", - " return -(4*(L - 2)**2 + 12*(L-2) + 8) / L**2\n", + " return -(4 * (L - 2) ** 2 + 12 * (L - 2) + 8) / L**2\n", + "\n", "\n", "L = 100\n", "print(f\"For L = {L}, We predict E = {E_prediction_all_the_same(L)}\")" @@ -223,15 +234,26 @@ " for k in range(M):\n", " # your code goes here\n", " pass\n", - " return E / (N*M)\n", - " \n", - "expected_values = [E_prediction_all_the_same(100), E_prediction_all_the_same(100), 0, '???'] \n", + " return E / (N * M)\n", "\n", - "f, axes = plt.subplots(ncols = 2, nrows = 2, figsize = (10,10))\n", + "\n", + "expected_values = [\n", + " E_prediction_all_the_same(100),\n", + " E_prediction_all_the_same(100),\n", + " 0,\n", + " \"???\",\n", + "]\n", + "\n", + "f, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))\n", "axes = axes.flatten()\n", "for ax, state, exp_value in zip(axes, states, expected_values):\n", - " show_state(state, ax = ax)\n", - " ax.text(0,-0.1, f\"Expected Value: {exp_value}, Your value: {energy(state)}\", transform = ax.transAxes)" + " show_state(state, ax=ax)\n", + " ax.text(\n", + " 0,\n", + " -0.1,\n", + " f\"Expected Value: {exp_value}, Your value: {energy(state)}\",\n", + " transform=ax.transAxes,\n", + " )" ] }, { @@ -254,6 +276,7 @@ "source": [ "## Solution\n", "\n", + "\n", "def energy(state):\n", " E = 0\n", " N, M = state.shape\n", @@ -262,22 +285,33 @@ " # handle the north and south neighbours\n", " for di in [1, -1]:\n", " if 0 <= (i + di) < N:\n", - " E -= state[i,j] * state[i+di, j]\n", - " \n", - " # handle the east and west neighbours\n", - " for dj in [1,-1]:\n", - " if 0 <= (j + dj) < M:\n", - " E -= state[i,j] * state[i, j+dj]\n", - " \n", - " return E / (N*M)\n", - " \n", - "expected_values = [E_prediction_all_the_same(100), E_prediction_all_the_same(100), 0, '???'] \n", + " E -= state[i, j] * state[i + di, j]\n", "\n", - "f, axes = plt.subplots(ncols = 2, nrows = 2, figsize = (10,10))\n", + " # handle the east and west neighbours\n", + " for dj in [1, -1]:\n", + " if 0 <= (j + dj) < M:\n", + " E -= state[i, j] * state[i, j + dj]\n", + "\n", + " return E / (N * M)\n", + "\n", + "\n", + "expected_values = [\n", + " E_prediction_all_the_same(100),\n", + " E_prediction_all_the_same(100),\n", + " 0,\n", + " \"???\",\n", + "]\n", + "\n", + "f, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))\n", "axes = axes.flatten()\n", "for ax, state, exp_value in zip(axes, states, expected_values):\n", - " show_state(state, ax = ax)\n", - " ax.text(0,-0.1, f\"Expected Value: {exp_value}, Your value: {energy(state)}\", transform = ax.transAxes)" + " show_state(state, ax=ax)\n", + " ax.text(\n", + " 0,\n", + " -0.1,\n", + " f\"Expected Value: {exp_value}, Your value: {energy(state)}\",\n", + " transform=ax.transAxes,\n", + " )" ] }, { @@ -315,9 +349,9 @@ } ], "source": [ - "L = 100 # How large the system should be\n", - "N = 100 # How many random samples to use\n", - "energies = [energy(np.random.choice([-1, 1], size = (L,L))) for _ in range(N)]\n", + "L = 100 # How large the system should be\n", + "N = 100 # How many random samples to use\n", + "energies = [energy(np.random.choice([-1, 1], size=(L, L))) for _ in range(N)]\n", "plt.hist(energies)\n", "print(f\"mean = {np.mean(energies)}, standard error = {np.std(energies) / np.sqrt(N)}\")" ] @@ -391,14 +425,21 @@ "def test_energy_function(energy_function):\n", " assert np.all(energy_function(state) == energy(state) for state in states)\n", "\n", + "\n", "def time_energy_function(energy_function):\n", " return [energy_function(state) for state in states]\n", "\n", + "\n", "def your_faster_energy_function(state):\n", - " return energy(state) # <-- replace this with your implementation and compare how fast it runs!\n", + " return energy(\n", + " state\n", + " ) # <-- replace this with your implementation and compare how fast it runs!\n", + "\n", "\n", "print(\"Naive baseline implementation\")\n", - "test_energy_function(energy) #this should always pass because it's just comparing to itself!\n", + "test_energy_function(\n", + " energy\n", + ") # this should always pass because it's just comparing to itself!\n", "naice = %timeit -o time_energy_function(energy)\n", "\n", "print(\"\\nYour version\")\n", @@ -443,16 +484,16 @@ } ], "source": [ - "\n", - "\n", "def numpy_energy(state):\n", - " E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) \n", - " return 2*E / np.product(state.shape)\n", + " E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])\n", + " return 2 * E / np.product(state.shape)\n", + "\n", + "\n", "test_energy_function(numpy_energy)\n", "\n", - "states[0][0,5] = -1\n", - "states[0][0,0] = -1\n", - "states[1][-1,-1] = 1\n", + "states[0][0, 5] = -1\n", + "states[0][0, 0] = -1\n", + "states[1][-1, -1] = 1\n", "print([energy(state) for state in states])\n", "print([energy2(state) for state in states])\n", "print([numba_energy(state) for state in states])\n", @@ -487,8 +528,10 @@ "def test_energy_function(energy_function):\n", " return [energy_function(state) for state in states]\n", "\n", + "\n", "from numba import jit\n", "\n", + "\n", "@jit(nopython=True)\n", "def numba_energy(state):\n", " E = 0\n", @@ -498,19 +541,20 @@ " # handle the north and south neighbours\n", " for di in [1, -1]:\n", " if 0 <= (i + di) < N:\n", - " E -= state[i,j] * state[i+di, j]\n", - " \n", + " E -= state[i, j] * state[i + di, j]\n", + "\n", " # handle the east and west neighbours\n", - " for dj in [1,-1]:\n", + " for dj in [1, -1]:\n", " if 0 <= (j + dj) < M:\n", - " E -= state[i,j] * state[i, j+dj]\n", - " \n", - " return E / (N*M)\n", + " E -= state[i, j] * state[i, j + dj]\n", + "\n", + " return E / (N * M)\n", + "\n", "\n", "def numpy_energy(state):\n", - " E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) \n", - " return 2*E / np.product(state.shape)\n", - " \n", + " E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])\n", + " return 2 * E / np.product(state.shape)\n", + "\n", "\n", "print(\"Naive baseline implementation\")\n", "naive = %timeit -o time_energy_function(energy)\n", @@ -573,16 +617,17 @@ " for j in range(M):\n", " # handle the north and south neighbours\n", " if 0 <= (i + 1) < N:\n", - " E -= state[i,j] * state[i+1, j]\n", - " \n", + " E -= state[i, j] * state[i + 1, j]\n", + "\n", " # handle the east and west neighbours\n", " if 0 <= (j + 1) < M:\n", - " E -= state[i,j] * state[i, j+1]\n", - " \n", - " return 2*E / (N*M)\n", + " E -= state[i, j] * state[i, j + 1]\n", + "\n", + " return 2 * E / (N * M)\n", + "\n", "\n", "print(\"\\nImproved Naive version\")\n", - "energy2(states[0]) #run the function once to let numba compile it before timing it\n", + "energy2(states[0]) # run the function once to let numba compile it before timing it\n", "e2 = %timeit -o test_energy_function(energy2)\n", "print(f\"Speedup: {naive.best/e2.best :.0f}x !\")" ] @@ -619,42 +664,45 @@ ], "source": [ "from numpy.random import default_rng\n", + "\n", "rng = default_rng(42)\n", "\n", - "def mcmc(initial_state, steps, T, energy = numba_energy):\n", + "\n", + "def mcmc(initial_state, steps, T, energy=numba_energy):\n", " N, M = initial_state.shape\n", " assert N == M\n", - " \n", + "\n", " current_state = initial_state.copy()\n", " E = N**2 * energy(state)\n", " for i in range(steps):\n", - " i, j = rng.integers(N, size = 2)\n", - " \n", + " i, j = rng.integers(N, size=2)\n", + "\n", " new_state = current_state.copy()\n", - " new_state[i,j] *= -1\n", + " new_state[i, j] *= -1\n", " new_E = N**2 * energy(new_state)\n", - " \n", + "\n", " if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform():\n", " current_state = new_state\n", " E = new_E\n", - " \n", - " return current_state \n", + "\n", + " return current_state\n", + "\n", "\n", "Ts = [4, 5, 50]\n", "\n", "ncols = 1 + len(Ts)\n", - "f, axes = plt.subplots(ncols = ncols, figsize = (5*ncols, 5))\n", + "f, axes = plt.subplots(ncols=ncols, figsize=(5 * ncols, 5))\n", "\n", - "show_state(initial_state, ax = axes[0])\n", - "axes[0].set(title = \"Initial state\")\n", + "show_state(initial_state, ax=axes[0])\n", + "axes[0].set(title=\"Initial state\")\n", "\n", "for T, ax in zip(Ts, axes[1:]):\n", " # initial_state = rng.choice([1,-1], size = (50,50))\n", - " initial_state = np.ones(shape = (50,50))\n", + " initial_state = np.ones(shape=(50, 50))\n", "\n", - " final_state = mcmc(initial_state, steps = 100_000, T = T) \n", - " show_state(final_state, ax = ax)\n", - " ax.set(title = f\"T = {T}\")" + " final_state = mcmc(initial_state, steps=100_000, T=T)\n", + " show_state(final_state, ax=ax)\n", + " ax.set(title=f\"T = {T}\")" ] }, { diff --git a/learning/02 Packaging it up.ipynb b/learning/02 Packaging it up.ipynb index fe94ff3..7721806 100644 --- a/learning/02 Packaging it up.ipynb +++ b/learning/02 Packaging it up.ipynb @@ -21,7 +21,9 @@ "\n", "# This loads some custom styles for matplotlib\n", "import json, matplotlib\n", - "with open(\"assets/matplotlibrc.json\") as f: matplotlib.rcParams.update(json.load(f))" + "\n", + "with open(\"assets/matplotlibrc.json\") as f:\n", + " matplotlib.rcParams.update(json.load(f))" ] }, { diff --git a/learning/assets/matplotlibrc.json b/learning/assets/matplotlibrc.json index e96b1b5..5609242 100644 --- a/learning/assets/matplotlibrc.json +++ b/learning/assets/matplotlibrc.json @@ -11,4 +11,4 @@ "svg.fonttype": "path", "figure.dpi": 100, "figure.figsize": [5,5] -} \ No newline at end of file +} diff --git a/requirements.txt b/requirements.txt index ae018e7..b8ccae8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,4 @@ ipykernel numpy scipy matplotlib -numba \ No newline at end of file +numba