Skip to content

Commit

Permalink
#24 Save and load checkpoint + misc changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@weka511
weka511 committed Jun 10, 2022
1 parent 4a7e59e commit d5898a9
Show file tree
Hide file tree
Showing 3 changed files with 583 additions and 464 deletions.
44 changes: 28 additions & 16 deletions book/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
from unittest import TestCase, main

class Geometry:
'''This class represents the apce in which the action occurs'''
'''This class represents the space in which the action occurs'''
def __init__(self,
L = array([1,1]),
sigma = 0.25,
Expand Down Expand Up @@ -50,8 +50,8 @@ def pack_disks(self, N=4):
def alloc(i,j):
'''Determine position of one disk'''
offset = 0 if i%2==0 else 0.5*Delta[1]
x0 = self.LowerBound[0] + (i+1) * Delta[0]
y0 = self.LowerBound[1] + offset + (j+1) * Delta[1]
x0 = self.LowerBound[0] + i * Delta[0]
y0 = self.LowerBound[1] + offset + j * Delta[1]
x1,y1 = self.move_to((x0,y0))
return (x1,y1)

Expand All @@ -69,21 +69,25 @@ def alloc(i,j):
coordinates = [alloc(i,j) for i in range(m) for j in range(n)]
return array(coordinates[0:N])

def create_Histograms(self,n=10):
return [Histogram(n) for _ in range(self.d)]
def create_Histograms(self,n=10,HistogramBins=zeros(0)):
self.HistogramBins = zeros((n,self.d)) if HistogramBins.size==0 else HistogramBins
return [Histogram(n,
h = self.HistogramBins[:,j]) for j in range(self.d)]

class Box(Geometry):
'''This class reprsents a simple box geometry without periodic boundary conditions'''
'''This class represents a simple box geometry without periodic boundary conditions'''
def __init__(self,
L = array([1,1]),
sigma = 0.25,
sigma = 0.125,
d = 2):
super().__init__(L = L, sigma=sigma, d = d)
super().__init__(L = L,
sigma = sigma,
d = d)
self.LowerBound = sigma*ones(d)
self.UpperBound = L - 2*sigma*ones(d)
self.UpperBound = L - sigma*ones(d)

def get_distance(self, X0,X1):
'''Calculate distance between two points using appropriate boundary conditions'''
'''Calculate Euclidean distance between two points'''
return norm(X0-X1)

def move_to(self,X):
Expand All @@ -96,14 +100,16 @@ class Torus(Geometry):
'''This class reprsents a box geometry with periodic boundary conditions, i.e. a torus.'''
def __init__(self,
L = array([1,1]),
sigma = 0.25,
sigma = 0.125,
d = 2):
super().__init__(L = L, sigma=sigma, d = d)
super().__init__(L = L,
sigma = sigma,
d = d)
self.LowerBound = zeros(d)
self.UpperBound = L

def get_distance(self, X0,X1):
'''Calculate distance between two points using appropriate boundary conditions'''
'''Calculate distance between two points using periodic boundary conditions'''
return norm(self.diff_vec(X0,X1))

def move_to(self,X):
Expand All @@ -124,7 +130,7 @@ def get_description(self):

def GeometryFactory(periodic = False,
L = array([1,1]),
sigma = 0.25,
sigma = 0.125,
d = 2):
'''Create a periodic or aperiodic Geometry'''
return Torus(L = L,
Expand All @@ -142,9 +148,10 @@ class Histogram:
def __init__(self,
n = 10,
x0 = 0,
xn = 1):
xn = 1,
h = zeros((0,0))):
self.n = n
self.h = [0] * n
self.h = zeros((n)) if h.size == 0 else h
self.x0 = x0
self.xn = xn

Expand All @@ -163,6 +170,11 @@ def get_hist(self):
bins = [self.x0 + i*(self.xn-self.x0)/self.n for i in range(self.n)]
return self.h,bins+[self.xn]

def bins(self):
for i in range(self.n):
yield(self.h[i])


class TestHistogram(TestCase):
def setUp(self):
self.histogram = Histogram()
Expand Down
91 changes: 64 additions & 27 deletions book/markov-disks.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,42 +18,40 @@
from argparse import ArgumentParser
from geometry import GeometryFactory
from matplotlib import rc
from matplotlib.pyplot import figure, legend, plot, savefig, show, title
from numpy import any, array, copy
from matplotlib.pyplot import bar, figure, legend, savefig, show, title
from numpy import any, array, copy, load, savez
from numpy.random import default_rng
from os.path import basename, splitext
from os.path import basename, exists, splitext
from shutil import copyfile

def markov_disks(X,
rng = default_rng(),
delta = array([0.01,0.01]),
geometry = GeometryFactory(),
sigma = 0.125):
geometry = GeometryFactory()):
'''Algorithm 2.9. Generating a hard disk configuration from an earlier valid configuration using MCMC'''

def can_move(k):
def can_move(k,X_proposed):
'''
Verify that proposed new position is within the geometry,
and that the resulting new configuration will be acceptable.
'''
if any(X[k,:] < geometry.LowerBound): return False
if any(geometry.UpperBound < X[k,:]): return False
if any(X_proposed < geometry.LowerBound): return False
if any(geometry.UpperBound < X_proposed): return False

for i in range(N):
if i!=k and geometry.get_distance(X[i,:],X[k,:])<2*sigma:
if i!=k and geometry.get_distance(X[i,:],X_proposed)<2*geometry.sigma:
return False

return True

N,d = X.shape
k = rng.integers(0,high=N)
Delta = -delta * 2* delta*rng.random(size=d)
x_save = copy(X[k,:]) # https://stackoverflow.com/questions/47181092/numpy-views-vs-copy-by-slicing
X[k,:] = geometry.move_to(X[k,:]+Delta) # Provisional move

if can_move(k):
X_proposed = geometry.move_to(X[k,:]+Delta)
if can_move(k,X_proposed):
X[k,:] = X_proposed
return k,X
else:
X[k,:] = x_save
return -1,X

def get_coordinate_description(coordinate):
Expand Down Expand Up @@ -113,27 +111,59 @@ def parse_arguments():
parser.add_argument('--frequency',
type = int,
default = 1000)
parser.add_argument('--restart',
action = 'store_true',
help = 'Restart from checkpoint')
return parser.parse_args()


class Checkpointer:
'''Used to save a configuration to a checkpoint, and restore from saved checkpoint'''
def __init__(self,file='check'):
self.path = f'{file}.npz'
self.backup = f'{self.path}~'

def load(self):
with load(self.path) as data:
X = data['X']
HistogramBins = data['HistogramBins']
return X,HistogramBins

def save(self,
X = [],
geometry = None):

if exists(self.path):
copyfile(self.path,self.backup)
savez(self.path,
X = X,
HistogramBins = geometry.HistogramBins)

if __name__=='__main__':
args = parse_arguments()
check = Checkpointer()
rng = default_rng()
delta = array(args.delta if len(args.delta)==args.d else args.delta * args.d)
geometry = GeometryFactory(periodic = args.periodic,
L = array(args.L if len(args.L)==args.d else args.L * args.d),
sigma = args.sigma,
d = args.d)
geometry = GeometryFactory(
periodic = args.periodic,
L = array(args.L if len(args.L)==args.d else args.L * args.d),
sigma = args.sigma,
d = args.d)
eta = geometry.get_density(N = args.Disks)
X = geometry.create_configuration(N=args.Disks)
n_accepted = 0
histograms = geometry.create_Histograms(n=args.bins)
if args.restart:
X,HistogramBins = check.load()
histograms = geometry.create_Histograms(n=args.bins,HistogramBins = HistogramBins)
else:
X = geometry.create_configuration(N = args.Disks)
histograms = geometry.create_Histograms(n=args.bins)

for epoch in range(args.N):
k,X = markov_disks(X,
rng = rng,
delta = delta,
geometry = geometry,
sigma = args.sigma)
k,X = markov_disks(
X,
rng = rng,
delta = delta,
geometry = geometry)

if epoch>args.burn:
if k>-1:
Expand All @@ -143,15 +173,22 @@ def parse_arguments():
histograms[i].add(x)
if epoch%args.frequency ==0:
print (f'Epoch {epoch:,}')
check.save(X = X,
geometry = geometry)


figure(figsize=(12,12))

for j in range(args.d):
h,bins = histograms[j].get_hist()
plot([0.5*(bins[i]+bins[i+1]) for i in range(len(h))],h,
label =f'{get_coordinate_description(j)}')
bar([0.5*(bins[i]+bins[i+1]) for i in range(len(h))],h,
width = [0.5*(bins[i+1]-bins[i]) for i in range(len(h))],
label = f'{get_coordinate_description(j)}',
alpha = 0.5,
color = 'xkcd:red' if j==0 else 'xkcd:blue' if j==1 else 'xkcd:green')
title(f'{geometry.get_description()} sigma = {args.sigma}, eta = {eta:.3f}, delta = {max(args.delta):.2g}, acceptance = {100*n_accepted/(args.N-args.burn):.3g}%')
legend()
savefig(get_plot_file_name(args.plot))

if args.show:
show()
Loading

0 comments on commit d5898a9

Please sign in to comment.