-
Notifications
You must be signed in to change notification settings - Fork 0
/
mcmc.py
57 lines (43 loc) · 1.67 KB
/
mcmc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import random
from math import log
import numpy as np
class IsingMCMC:
"""Ising MCMC with simulated annealing."""
def __init__(self, graph, interaction, fields, temperature):
self.J = interaction # interaction potential
self.h = fields # local fields
self.B = 1 / temperature
self.graph = graph
def run(self, init_states, steps, gamma=.01):
"""Run the simulation.
Args:
init_states: Dictionary mapping nodes to their initial state.
steps: Number of steps.
temperature: Initial temperature.
gamma: Cooling factor.
Returns:
The final states dictionary and the list of energies by time.
"""
states = init_states.copy()
beta = self.B
energies = np.zeros(steps + 1)
energies[0] = self.energy(states)
nodes = list(self.graph.nodes)
for t in range(steps):
for node in nodes:
state = states[node]
nstates = [states[n] for n in self.graph.neighbors(node)]
# Metropolis–Hastings
delta = 2 * self.J * state * np.sum(nstates) + 2 * self.h[node] * state
if random.random() < np.exp(- beta * delta):
states[node] *= -1 # flip the spin
random.shuffle(nodes)
beta *= 1 + log(1 + gamma * t) # simulated annealing
energies[1 + t] = self.energy(states)
return states, energies
def energy(self, states):
"""Calculate the system energy."""
E = 0
for (u, v) in self.graph.edges:
E -= self.J * states[u] * states[v]
return E