globally remove the 1/L**2 factor from definitions

code/src/MCFF/ising_model.py

@@ -75,7 +75,7 @@ def energy(state):
             if 0 <= (j + 1) < M:
                 E -= state[i, j] * state[i, j + 1]
 
-    return 2 * E / (N * M)
+    return 2 * E
 
 
 def energy_numpy(state):
@@ -93,7 +93,7 @@ def energy_numpy(state):
         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)
+    return 2 * E
 
 
 @jit(nopython=True, nogil=True)
@@ -108,4 +108,4 @@ def energy_difference(state, site):
             0 <= (j + dj) < M
         ):  # ignore neighbours not in the NxN grid
             h += state[i + di, j + dj]
-    return 4 * state[i, j] * h / (N * M)
+    return 4 * state[i, j] * h

code/src/MCFF/mcmc.py

@@ -7,9 +7,11 @@ from numba import jit
 from .ising_model import energy, energy_difference
 
 
-def mcmc(initial_state, steps, T, rng=np.random.default_rng()):
+def mcmc_original(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`.
 
+    This is the original version of the algorithm, written in notebook 03.
+
     Parameters
     ----------
     initial_state : array_like
@@ -30,13 +32,13 @@ def mcmc(initial_state, steps, T, rng=np.random.default_rng()):
     assert N == M
 
     current_state = initial_state.copy()
-    E = N**2 * energy(initial_state)
+    E = energy(initial_state)
     for i in range(steps):
         i, j = rng.integers(N, size=2)
 
         new_state = current_state.copy()
         new_state[i, j] *= -1
-        new_E = N**2 * energy(new_state)
+        new_E = energy(new_state)
 
         if (new_E < E) or np.exp(-(new_E - E) / T) > rng.uniform():
             current_state = new_state
@@ -49,6 +51,27 @@ def mcmc(initial_state, steps, T, rng=np.random.default_rng()):
 def mcmc_generator(
     initial_state, steps, T, stepsize=1000, energy_difference=energy_difference
 ):
+    r"""Yields the state after each MCMC step. Performs `steps` monte carlo steps each of size `stepsize` starting from `initial_state` at temperature `T`.
+
+    Parameters
+    ----------
+    initial_state : array_like
+        Numpy array of shape (N,N) with values +1,-1
+    steps : int
+        The number of steps to take.
+    T : float
+        The temperature at which to perform the simulation.
+    stepsize : int, default = 1000
+        How many times to try to flip a pixel between for each step.
+
+    energy_difference : function(state, site), default = ising_model.energy_difference
+        A function that gives the energy change from flippings `site = (i, j)`of `state`.
+
+    Returns
+    -------
+    np.array
+        The final state.
+    """
     N, M = initial_state.shape
     assert N == M
 
@@ -58,7 +81,7 @@ def mcmc_generator(
             i, j = np.random.randint(N), np.random.randint(N)
 
             # calculate the energy change if we were to flip this pixel but don't actually do it
-            change_in_E = N**2 * energy_difference(current_state, (i, j))
+            change_in_E = energy_difference(current_state, (i, j))
 
             if change_in_E < 0 or np.exp(-change_in_E / T) > np.random.random():
                 current_state[i, j] *= -1  # accept the change!

code/tests/test_energy.py

@@ -18,7 +18,7 @@ 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)
 
 
 def test_exact_energies():

learning/04 Testing.ipynb

@@ -64,7 +64,7 @@
     "# Note that only functions whose name begins with test_ get run by pytest\n",
     "def E_prediction_all_the_same(L): \n",
     "    \"The exact energy in for the case where all spins are up or down\"\n",
-    "    return -(4*(L - 2)**2 + 12*(L-2) + 8) / L**2\n",
+    "    return -(4*(L - 2)**2 + 12*(L-2) + 8)\n",
     "\n",
     "def test_exact_energies():\n",
     "    for state in [all_up, all_down]:\n",

learning/05 Adding functionality.ipynb

@@ -65,7 +65,7 @@
     "    ...\n",
     "\n",
     "    current_state = initial_state.copy()\n",
-    "    E = N**2 * energy(current_state)\n",
+    "    E = energy(current_state)\n",
     "    for i in range(steps):\n",
     "        measurements[i] = measurement(state)\n",
     "        ...\n",
@@ -127,7 +127,7 @@
     {
      "data": {
       "text/plain": [
-       "(<generator object my_range at 0x7fab3acae4a0>, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
+       "(<generator object my_range at 0x7fe3b2359660>, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
       ]
      },
      "execution_count": 2,
@@ -158,7 +158,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "id": "f73d6335-6514-45b1-9128-d72122d8b0b7",
    "metadata": {},
    "outputs": [
@@ -166,15 +166,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "19.3 ms ± 156 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "18.4 ms ± 364 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
-      "10.8 ms ± 869 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "9.51 ms ± 309 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
       "0x slowdown!\n"
      ]
     }
    ],
    "source": [
     "from MCFF.ising_model import energy, all_up_state\n",
-    "from MCFF.mcmc import mcmc\n",
+    "from MCFF.mcmc import mcmc_original\n",
     "\n",
     "\n",
     "@jit(nopython=True, nogil=True)\n",
@@ -183,14 +183,14 @@
     "    assert N == M\n",
     "\n",
     "    current_state = initial_state.copy()\n",
-    "    E = N**2 * energy(current_state)\n",
+    "    E = energy(current_state)\n",
     "    for _ in range(steps):\n",
     "        for _ in range(stepsize):\n",
     "            i, j = np.random.randint(N), np.random.randint(N)\n",
     "\n",
     "            # modify the state a little, here we just flip a random pixel\n",
     "            current_state[i, j] *= -1\n",
-    "            new_E = N**2 * energy(current_state)\n",
+    "            new_E = energy(current_state)\n",
     "\n",
     "            if (new_E < E) or np.exp(-(new_E - E) / T) > np.random.random():\n",
     "                E = new_E\n",
@@ -203,7 +203,7 @@
     "N_steps = 1000\n",
     "stepsize = 1\n",
     "initial_state = all_up_state(20)\n",
-    "without_yield = %timeit -o mcmc(initial_state, steps = N_steps, T = 5)\n",
+    "without_yield = %timeit -o mcmc_original(initial_state, steps = N_steps, T = 5)\n",
     "with_yield = %timeit -o [np.mean(s) for s in mcmc_generator(initial_state, T = 5, steps = N_steps, stepsize = 1)]\n",
     "print(f\"{with_yield.best / without_yield.best:.0f}x slowdown!\")"
    ]