1
\$\begingroup\$

I've written this code to simulate the Ising Model. The idea is a system represented as a matrix with each element having value 1 or -1. The input is an initial matrix of 1 and -1; it is then being subjected to this algorithm:

  1. Pick a random element of the matrix
  2. calculate the energy change of the system if that element is flipped (turn 1 to -1 vice versa)
  3. if the flipping reduces total energy (after/before < 0): do it, else: flip it with probability

$$ P = \exp(-\text{constant}*\text{energy difference}) $$ $$0 \le P \le 1$$

(I added an extra term in the code to modify its dynamics.)

Energy is calculated by observing the spin of an element with its neighbor, so I use two matrices for this calculation:

  1. the main matrix (the one subjected to the function),
  2. an adjacency matrix (representing connection between each elements from the main matrix: 1 if connected, else 0).

The idea is to multiply the row from the adjacency matrix to all elements of the main matrix.

The code works but since I have to apply the function for thousands of times and it uses two (rather large) matrices, the calculation takes a lot of time. Perhaps ideas for speed-ups (or just for writing a better code in general) and possibility to use numba which currently gives me complicated error messages.

import numpy as np

n = 64
n2 = n**2

main_matrix = 2*np.random.randint(2, size=(n, n))-1

adj_matrix = np.random.randint(2, size=(n2, n2)) # Assume random network is used

def mcmove(beta):
  a = np.random.randint(0, n)
  b = np.random.randint(0, n)
  s = main_matrix[a, b]

  energy_difference = np.sum(np.multiply(adj_matrix[a*n + b], main_matrix.reshape(1,-1)[0]))/np.sum(adj_matrix[a*n + b])
  cost = 2*s*(energy_difference - s*abs(np.sum(main_matrix))/n2)

  if cost < 0:
    s *= -1
  elif np.random.rand() < np.exp(-cost*beta):
    s *= -1
  main_matrix[a, b] = s
  return main_matrix

# How this function is used
total_spin = []

for i in range(1000):
  for i in range(n2):
    mcmove(0.001*i)
  total_spin.append(np.sum(main_matrix)) # sum all spins

total_spin
\$\endgroup\$

1 Answer 1

1
\$\begingroup\$

Stop using np.random.randint and family; they're deprecated in favour of the new generators.

You can assign a and b in one go by generating an integer array of length 2.

Add PEP484 type hints.

Don't np.multiply; in this case use np.dot since effectively you are calculating a dot product.

Don't reshape; your reshape is equivalent to a flatten().

Use np.abs instead of abs.

Don't elif; this is equivalent to an or sub-expression.

Don't return main_matrix; that's ignored and you're mutating in-place.

Use a different letter for your inner for i iteration, probably j.

Rather than 0.001 * i, find this value directly by iterating over a call to np.linspace.

Rather than total_spin.append, pre-allocate a properly typed, empty numpy array.

Other than the above there isn't much you can improve in terms of performance since the iterations are not independent: an evolved main_matrix depends on the value of main_matrix from the previous iteration. At that point you would drop down to a better language like C and calls to cblas or equivalent.

Suggested

import numpy as np

n = 64
n2 = n**2

rand = np.random.default_rng(seed=0)
main_matrix = rand.choice((-1, 1), size=(n, n))
adj_matrix = rand.integers(2, size=(n2, n2))  # Assume random network is used


def mcmove(beta: float) -> None:
    a, b = rand.integers(low=0, high=n, size=2)
    s = main_matrix[a, b]

    adj_row = adj_matrix[a*n + b, :]
    energy_difference = np.dot(adj_row, main_matrix.flatten()) / adj_row.sum()
    cost = 2 * s * (energy_difference - s * np.abs(main_matrix.sum()) / n2)

    if cost < 0 or rand.random() < np.exp(-cost * beta):
        main_matrix[a, b] = -s


def main() -> None:
    # How this function is used
    N = 1000
    total_spin = np.empty(shape=N, dtype=int)

    for i in range(N):
        for j in np.linspace(0, (n2 - 1)/N, n2):
            mcmove(j)
        total_spin[i] = main_matrix.sum()  # sum all spins

    print(total_spin)


if __name__ == '__main__':
    main()
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.