tom/ReCoDE_MCMCFF
1
0
Fork 0
mirror of https://github.com/ImperialCollegeLondon/ReCoDE_MCMCFF.git synced 2025-06-26 08:51:16 +02:00
Code Issues Packages Projects Releases Wiki Activity

run black on everything

Browse Source 
autoformatting
This commit is contained in:
Tom Hodson 2022-06-09 09:57:58 +02:00
parent 9de5c8c303
commit 9eb933cb14
12 changed files with 183 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
code/src/MCFF/ising_model.py
View File

@ -7,15 +7,18 @@ States are represented by NxM numpy arrays with values +/-1. These functions are
import numpy as np import numpy as np
from numba import jit from numba import jit
def all_up_state(N): def all_up_state(N):
"The NxN all up state of the Ising model." "The NxN all up state of the Ising model."
return np.ones([N, N]) return np.ones([N, N])
def all_down_state(N): def all_down_state(N):
"The NxN all down state of the Ising model." "The NxN all down state of the Ising model."
return -np.ones([N, N]) 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. r"""Return the a random NxN state of the Ising model.
Parameters Parameters
@ -32,13 +35,13 @@ def random_state(N, magnetisation = 0.5, rng = np.random.default_rng()):
np.array np.array
A random NxN state. 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 # 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. # 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 # 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): def energy(state):
r"""Compute the energy of a state of the Ising Model with open boundary conditions. r"""Compute the energy of a state of the Ising Model with open boundary conditions.
@ -57,13 +60,14 @@ def energy(state):
for j in range(M): for j in range(M):
# handle the north and south neighbours # handle the north and south neighbours
if 0 <= (i + 1) < N: 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 # handle the east and west neighbours
if 0 <= (j + 1) < M: if 0 <= (j + 1) < M:
E -= state[i,j] * state[i, j+1] E -= state[i, j] * state[i, j + 1]
return 2 * E / (N * M)
return 2*E / (N*M)
def energy_numpy(state): def energy_numpy(state):
r"""Compute the energy of a state of the Ising Model with open boundary conditions. r"""Compute the energy of a state of the Ising Model with open boundary conditions.
@ -79,6 +83,5 @@ def energy_numpy(state):
float64 float64
The interaction energy per site. The interaction energy per site.
""" """
E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])
return 2*E / np.product(state.shape) return 2 * E / np.product(state.shape)

7
code/src/MCFF/mcmc.py
View File

@ -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`. r"""Return the state after performing `steps` monte carlo steps starting from `initial_state` at temperature `T`.
Parameters Parameters
@ -28,10 +29,10 @@ def mcmc(initial_state, steps, T, rng = np.random.default_rng()):
current_state = initial_state.copy() current_state = initial_state.copy()
E = N**2 * energy(state) E = N**2 * energy(state)
for i in range(steps): 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 = current_state.copy()
new_state[i,j] *= -1 new_state[i, j] *= -1
new_E = N**2 * energy(new_state) new_E = N**2 * energy(new_state)
if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform(): if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform():

15
code/tests/test_energy.py
View File

@ -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 from MCFF.ising_model import energy, energy_numpy
# create some test states # create some test states
all_up = np.ones([100,100]) all_up = np.ones([100, 100])
all_down = -np.ones([100,100]) all_down = -np.ones([100, 100])
random = np.random.choice([-1, 1], size = (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 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] states = [all_up, all_down, random, custom]
def E_prediction_all_the_same(L): def E_prediction_all_the_same(L):
"The exact energy in for the case where all spins are up or down" "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(): def test_exact_energies():
for state in [all_up, all_down]: for state in [all_up, all_down]:
L = state.shape[0] L = state.shape[0]
assert energy(state) == E_prediction_all_the_same(L) assert energy(state) == E_prediction_all_the_same(L)
def test_energy_implementations(): def test_energy_implementations():
for state in states: for state in states:
assert np.allclose(energy(state), energy_numpy(state)) assert np.allclose(energy(state), energy_numpy(state))

11
code/tests/test_energy_using_hypothesis.py
View File

@ -5,8 +5,13 @@ from hypothesis.extra import numpy as hnp
from MCFF.ising_model import energy, energy_numpy from MCFF.ising_model import energy, energy_numpy
@given(hnp.arrays(dtype = int,
shape = hnp.array_shapes(min_dims = 2, max_dims = 2), @given(
elements = st.sampled_from([1, -1]))) 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): def test_generated_states(state):
assert np.allclose(energy(state), energy_numpy(state)) assert np.allclose(energy(state), energy_numpy(state))

200
learning/01 Introduction.ipynb
View File

@ -51,9 +51,13 @@
"\n", "\n",
"# This loads some custom styles for matplotlib\n", "# This loads some custom styles for matplotlib\n",
"import json, matplotlib\n", "import json, matplotlib\n",
"with open(\"assets/matplotlibrc.json\") as f: matplotlib.rcParams.update(json.load(f))\n",
"\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": [ "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", "\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)" "show_state(state)"
] ]
}, },
@ -132,18 +139,21 @@
} }
], ],
"source": [ "source": [
"all_up = np.ones([100,100])\n", "all_up = np.ones([100, 100])\n",
"all_down = -np.ones([100,100])\n", "all_down = -np.ones([100, 100])\n",
"random = np.random.choice([-1, 1], size = (100,100))\n", "random = np.random.choice([-1, 1], size=(100, 100))\n",
"\n", "\n",
"from matplotlib.image import imread\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", "\n",
"states = [all_up, all_down, random, custom]\n", "states = [all_up, all_down, random, custom]\n",
"\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", "for ax, state in zip(axes, states):\n",
" show_state(state, ax = ax)" " show_state(state, ax=ax)"
] ]
}, },
{ {
@ -181,7 +191,8 @@
], ],
"source": [ "source": [
"def E_prediction_all_the_same(L):\n", "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", "\n",
"L = 100\n", "L = 100\n",
"print(f\"For L = {L}, We predict E = {E_prediction_all_the_same(L)}\")" "print(f\"For L = {L}, We predict E = {E_prediction_all_the_same(L)}\")"
@ -223,15 +234,26 @@
" for k in range(M):\n", " for k in range(M):\n",
" # your code goes here\n", " # your code goes here\n",
" pass\n", " pass\n",
" return E / (N*M)\n", " return E / (N * M)\n",
" \n",
"expected_values = [E_prediction_all_the_same(100), E_prediction_all_the_same(100), 0, '???'] \n",
"\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", "axes = axes.flatten()\n",
"for ax, state, exp_value in zip(axes, states, expected_values):\n", "for ax, state, exp_value in zip(axes, states, expected_values):\n",
" show_state(state, ax = ax)\n", " show_state(state, ax=ax)\n",
" ax.text(0,-0.1, f\"Expected Value: {exp_value}, Your value: {energy(state)}\", transform = ax.transAxes)" " 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": [ "source": [
"## Solution\n", "## Solution\n",
"\n", "\n",
"\n",
"def energy(state):\n", "def energy(state):\n",
" E = 0\n", " E = 0\n",
" N, M = state.shape\n", " N, M = state.shape\n",
@ -262,22 +285,33 @@
" # handle the north and south neighbours\n", " # handle the north and south neighbours\n",
" for di in [1, -1]:\n", " for di in [1, -1]:\n",
" if 0 <= (i + di) < N:\n", " if 0 <= (i + di) < N:\n",
" E -= state[i,j] * state[i+di, j]\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",
"\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", "axes = axes.flatten()\n",
"for ax, state, exp_value in zip(axes, states, expected_values):\n", "for ax, state, exp_value in zip(axes, states, expected_values):\n",
" show_state(state, ax = ax)\n", " show_state(state, ax=ax)\n",
" ax.text(0,-0.1, f\"Expected Value: {exp_value}, Your value: {energy(state)}\", transform = ax.transAxes)" " 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": [ "source": [
"L = 100 # How large the system should be\n", "L = 100 # How large the system should be\n",
"N = 100 # How many random samples to use\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", "energies = [energy(np.random.choice([-1, 1], size=(L, L))) for _ in range(N)]\n",
"plt.hist(energies)\n", "plt.hist(energies)\n",
"print(f\"mean = {np.mean(energies)}, standard error = {np.std(energies) / np.sqrt(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", "def test_energy_function(energy_function):\n",
" assert np.all(energy_function(state) == energy(state) for state in states)\n", " assert np.all(energy_function(state) == energy(state) for state in states)\n",
"\n", "\n",
"\n",
"def time_energy_function(energy_function):\n", "def time_energy_function(energy_function):\n",
" return [energy_function(state) for state in states]\n", " return [energy_function(state) for state in states]\n",
"\n", "\n",
"\n",
"def your_faster_energy_function(state):\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", "\n",
"print(\"Naive baseline implementation\")\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", "naice = %timeit -o time_energy_function(energy)\n",
"\n", "\n",
"print(\"\\nYour version\")\n", "print(\"\\nYour version\")\n",
@ -443,16 +484,16 @@
} }
], ],
"source": [ "source": [
"\n",
"\n",
"def numpy_energy(state):\n", "def numpy_energy(state):\n",
" E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) \n", " E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])\n",
" return 2*E / np.product(state.shape)\n", " return 2 * E / np.product(state.shape)\n",
"\n",
"\n",
"test_energy_function(numpy_energy)\n", "test_energy_function(numpy_energy)\n",
"\n", "\n",
"states[0][0,5] = -1\n", "states[0][0, 5] = -1\n",
"states[0][0,0] = -1\n", "states[0][0, 0] = -1\n",
"states[1][-1,-1] = 1\n", "states[1][-1, -1] = 1\n",
"print([energy(state) for state in states])\n", "print([energy(state) for state in states])\n",
"print([energy2(state) for state in states])\n", "print([energy2(state) for state in states])\n",
"print([numba_energy(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", "def test_energy_function(energy_function):\n",
" return [energy_function(state) for state in states]\n", " return [energy_function(state) for state in states]\n",
"\n", "\n",
"\n",
"from numba import jit\n", "from numba import jit\n",
"\n", "\n",
"\n",
"@jit(nopython=True)\n", "@jit(nopython=True)\n",
"def numba_energy(state):\n", "def numba_energy(state):\n",
" E = 0\n", " E = 0\n",
@ -498,19 +541,20 @@
" # handle the north and south neighbours\n", " # handle the north and south neighbours\n",
" for di in [1, -1]:\n", " for di in [1, -1]:\n",
" if 0 <= (i + di) < N:\n", " if 0 <= (i + di) < N:\n",
" E -= state[i,j] * state[i+di, j]\n", " E -= state[i, j] * state[i + di, j]\n",
" \n", "\n",
" # handle the east and west neighbours\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", " if 0 <= (j + dj) < M:\n",
" E -= state[i,j] * state[i, j+dj]\n", " E -= state[i, j] * state[i, j + dj]\n",
" \n", "\n",
" return E / (N*M)\n", " return E / (N * M)\n",
"\n",
"\n", "\n",
"def numpy_energy(state):\n", "def numpy_energy(state):\n",
" E = - np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:]) \n", " E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])\n",
" return 2*E / np.product(state.shape)\n", " return 2 * E / np.product(state.shape)\n",
" \n", "\n",
"\n", "\n",
"print(\"Naive baseline implementation\")\n", "print(\"Naive baseline implementation\")\n",
"naive = %timeit -o time_energy_function(energy)\n", "naive = %timeit -o time_energy_function(energy)\n",
@ -573,16 +617,17 @@
" for j in range(M):\n", " for j in range(M):\n",
" # handle the north and south neighbours\n", " # handle the north and south neighbours\n",
" if 0 <= (i + 1) < N:\n", " if 0 <= (i + 1) < N:\n",
" E -= state[i,j] * state[i+1, j]\n", " E -= state[i, j] * state[i + 1, j]\n",
" \n", "\n",
" # handle the east and west neighbours\n", " # handle the east and west neighbours\n",
" if 0 <= (j + 1) < M:\n", " if 0 <= (j + 1) < M:\n",
" E -= state[i,j] * state[i, j+1]\n", " E -= state[i, j] * state[i, j + 1]\n",
" \n", "\n",
" return 2*E / (N*M)\n", " return 2 * E / (N * M)\n",
"\n",
"\n", "\n",
"print(\"\\nImproved Naive version\")\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", "e2 = %timeit -o test_energy_function(energy2)\n",
"print(f\"Speedup: {naive.best/e2.best :.0f}x !\")" "print(f\"Speedup: {naive.best/e2.best :.0f}x !\")"
] ]
@ -619,42 +664,45 @@
], ],
"source": [ "source": [
"from numpy.random import default_rng\n", "from numpy.random import default_rng\n",
"\n",
"rng = default_rng(42)\n", "rng = default_rng(42)\n",
"\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", " N, M = initial_state.shape\n",
" assert N == M\n", " assert N == M\n",
" \n", "\n",
" current_state = initial_state.copy()\n", " current_state = initial_state.copy()\n",
" E = N**2 * energy(state)\n", " E = N**2 * energy(state)\n",
" for i in range(steps):\n", " for i in range(steps):\n",
" i, j = rng.integers(N, size = 2)\n", " i, j = rng.integers(N, size=2)\n",
" \n", "\n",
" new_state = current_state.copy()\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", " new_E = N**2 * energy(new_state)\n",
" \n", "\n",
" if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform():\n", " if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform():\n",
" current_state = new_state\n", " current_state = new_state\n",
" E = new_E\n", " E = new_E\n",
" \n", "\n",
" return current_state \n", " return current_state\n",
"\n",
"\n", "\n",
"Ts = [4, 5, 50]\n", "Ts = [4, 5, 50]\n",
"\n", "\n",
"ncols = 1 + len(Ts)\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", "\n",
"show_state(initial_state, ax = axes[0])\n", "show_state(initial_state, ax=axes[0])\n",
"axes[0].set(title = \"Initial state\")\n", "axes[0].set(title=\"Initial state\")\n",
"\n", "\n",
"for T, ax in zip(Ts, axes[1:]):\n", "for T, ax in zip(Ts, axes[1:]):\n",
" # initial_state = rng.choice([1,-1], size = (50,50))\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", "\n",
" final_state = mcmc(initial_state, steps = 100_000, T = T) \n", " final_state = mcmc(initial_state, steps=100_000, T=T)\n",
" show_state(final_state, ax = ax)\n", " show_state(final_state, ax=ax)\n",
" ax.set(title = f\"T = {T}\")" " ax.set(title=f\"T = {T}\")"
] ]
}, },
{ {

4
learning/02 Packaging it up.ipynb
View File

@ -21,7 +21,9 @@
"\n", "\n",
"# This loads some custom styles for matplotlib\n", "# This loads some custom styles for matplotlib\n",
"import json, 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))"
] ]
}, },
{ {