import numpy as np
import scipy.stats as scs
import matplotlib.pyplot as plt
import numpy.random as rand





########################
# Ising configurations #
########################

class IsingConfiguration:
	"""
	A container for configurations of spins +/- 1 on a grid. IsingConfiguration(N, mat) returns the configuration of size n and with matrix mat. If no parameter mat is specified, the spins of the configuration are random.

	Useful methods:
	- plot()
	- markov_plot(beta)
	- energy(): computes the sum over neighbors of products of spins.
	- magnetization(): computes the sum of spins of the configuration.
	- random_sites(): picks random sites on the grid.
	- switch_sites(y, u, beta): switches the sites y according to the parameters u and beta.
	- neighbors(i, j): the list of neighbors of the site with coordinates (i,j).
	- transform(n, beta): performs n transitions of the Markov chain with parameter beta.
	"""
	def __init__(self, N, mat=None):
		self.size = N
		if type(mat) == np.ndarray:
			self.matrix = mat
		else:
			self.matrix = 2*scs.bernoulli(0.5).rvs(size=(self.size, self.size)) - 1

	def __repr__(self):
		return f"Ising configuration on a {self.size}*{self.size} grid"

	def __call__(self, i, j):
		return self.matrix[i,j]

	def __eq__(self, other):
		return ((self.size == other.size) and (self.matrix == other.matrix))

	def plot(self):
		fig, ax = plt.subplots(figsize=(10,10))
		for i in range(self.size):
			for j in range(self.size):
				if self.matrix[i,j] == 1:
					ax.fill([i, i+1, i+1, i], [j, j, j+1, j+1], color=(0.9,0.3,0.3))
				else:
					ax.fill([i, i+1, i+1, i], [j, j, j+1, j+1], color=(0.7,0.7,0.7))
		ax.set_axis_off()
		ax.set_aspect(1)
		plt.show()

	def energy(self):
		"""
		The total energy of the system.
		"""
		return -int(sum([self(i,j) * self(i+1,j) for i in range(self.size-1) for j in range(self.size)]) + sum([self(i,j) * self(i,j+1) for i in range(self.size) for j in range(self.size-1)]))
		
	def magnetization(self):
		"""
		The sum of all spins.
		"""
		return int(np.sum(self.matrix))

	def empty(self):
		"""
		Sets all spins to -1.
		"""
		self.matrix = - np.ones([self.size, self.size], dtype=int)

	def fill(self):
		"""
		Sets all spins to +1.
		"""
		self.matrix = + np.ones([self.size, self.size], dtype=int)

	def fix_boundary(self, boundary="positive"):
		"""
		Fix the spins on the boundary.
		"""
		if boundary=="positive":
			for i in range(self.size):
				self.matrix[0, i] = +1
				self.matrix[self.size-1, i] = +1
				self.matrix[i, 0] = +1
				self.matrix[i, self.size-1] = +1
		elif boundary=="war":
			for i in range(self.size):
				self.matrix[0, i] = +1
				self.matrix[self.size-1, i] = -1
		else:
			pass
		
	def random_sites(self, boundary=None, size=1):
		"""
		Picks sites at random (with an optional boundary condition, which
		forbids a certain region of the grid).
		"""

		if boundary=="positive":	
			return  rand.randint(1, self.size-1, size=[size,2])
		elif boundary=="war":	
			return np.array([rand.randint(1, self.size-1, size), rand.randint(0, self.size, size)]).transpose()
		else:
			return rand.randint(0, self.size, size=[size,2])
			
	def neighbors(self, i, j):
		"""
		Returns the list of neighbors of (i,j).
		"""
		res = []
		if i < self.size-1:
			res.append((i+1, j))
		if j < self.size-1:
			res.append((i, j+1))
		if i > 0:
			res.append((i-1, j))
		if j > 0:
			res.append((i, j-1))
		return res

	def switch_site(self, y, u, beta=0):
		"""
		Upgrade the site y by comparing u to a threshold, which is a function of beta.
		"""
		num = np.exp(2*beta*sum([self.matrix[x[0],x[1]] for x in self.neighbors(y[0], y[1])])) 
		p = num/(1+num)
		if u < p:
			self.matrix[y[0], y[1]] = +1
		else:
			self.matrix[y[0], y[1]] = -1

	def switch_sites(self, y, u, beta=0, reverse=False):
		"""
		Upgrade all the sites in y by comparing the values in u to certain thresholds computed in terms of beta.
		"""
		if reverse:
			for i in range(len(y)-1, -1, -1):
				self.switch_site(y[i], u[i], beta)
		else:
			for i in range(len(y)):
				self.switch_site(y[i], u[i], beta)

	def transform(self, n=1, beta=0):
		"""
		Performs n transition of the Markov chain with parameter beta.
		"""
		y = self.random_sites(size=n)
		u = rand.random(size=n)
		self.switch_sites(y, u, beta)


def MarkovIsingConfiguration(N, n, beta=0, boundary=None):
	"""
	Computes a random Ising configuration with size N and parameter beta, by performing n transitions of a Markov chain.

	Parameters:
	- N: the size of the grid.
	- n: the number of transitions of the Markov chain.
	- beta: the parameter of the probability measure.
	- boundary: whether the boundary is fixed ("positive" or "war") or not (None).
	"""
	I = IsingConfiguration(N)
	I.fix_boundary(boundary)
	y = I.random_sites(boundary, size=n)
	u = rand.random(size=n)
	I.switch_sites(y, u, beta)
	return I				


def EmptyIsingConfiguration(N, boundary=None):
	res = IsingConfiguration(N, mat=-np.ones([N,N], dtype=int))
	res.fix_boundary(boundary)
	return res


def FullIsingConfiguration(N, boundary=None):
	res = IsingConfiguration(N, mat=np.ones([N,N], dtype=int))
	res.fix_boundary(boundary)
	return res


def RandomIsingConfiguration(N, beta=0, boundary=None):
	"""
	Computes a random Ising configuration with size N and parameter beta, by using the coupling from the past method.

	Parameters:
	- N: the size of the grid.
	- beta: the parameter of the probability measure.
	- boundary: whether the boundary is fixed ("positive" or "war") or not (None).
	"""
	count = 0
	L = []
	M = []
	EI = EmptyIsingConfiguration(N, boundary)
	FI = FullIsingConfiguration(N, boundary)
	test = True
	while test:
		L += [EI.random_sites(boundary)[0] for k in range(2**count)]
		M += [rand.random() for k in range(2**count)]
		count += 1
		EI.switch_sites(L, M, beta, True)
		FI.switch_sites(L, M, beta, True)
		if EI.magnetization() == FI.magnetization():
			test = False
		else:
			EI = EmptyIsingConfiguration(N, boundary)
			FI = FullIsingConfiguration(N, boundary)
	return EI













####################
# Plane partitions #
####################



class HexaPoint:
	"""
	A container for the vertices of the hexagons of an hexagonal tiling (honeycomb).

	Useful methods:
	- graphic_coordinates()
	- draw_on_ax(ax)
	- neighbors(N): compute the list of neighbors in the honeycomb of size N.
	"""
	def __init__(self, a, b):
		self.x = a
		self.y = b

	def __eq__(self, other):
		return ((self.x == other.x) and (self.y == other.y))

	def graphic_coordinates(self, rotate=True):
		a = self.x
		b = self.y
		if rotate:
			if (b >= 0) and (b % 2 == 0):
				return [np.sqrt(3)*a/2, (a + b)/2]
			elif (b > 0) and (b % 2 == 1):
				return [np.sqrt(3)*a/2 + 1/(2*np.sqrt(3)), (a + b)/2]
			elif (b < 0) and (b % 2 == 1):
				return [np.sqrt(3)*(a - b/2)/2 - 1/(4*np.sqrt(3)), (a + b/2)/2 - 1/4]
			elif (b < 0) and (b % 2 == 0):
				return [np.sqrt(3)*(a - b/2)/2, (a + b/2)/2]
		else:
			if (b >= 0) and (b % 2 == 0):
				return [a + b/4, np.sqrt(3)*b/4]
			elif (b > 0) and (b % 2 == 1):
				return [a + b/4 + 1/4, np.sqrt(3)*b/4 - 1/(4*np.sqrt(3))]
			elif (b < 0) and (b % 2 == 1):
				return [a - b/4 - 1/4, np.sqrt(3)*b/4 - 1/(4*np.sqrt(3))]
			elif (b < 0) and (b % 2 == 0):
				return [a - b/4, np.sqrt(3)*b/4]

	def neighbors(self, N):
		res = []
		a = self.x
		b = self.y
		if b in np.arange(0, 2*N, 2):
			if 0 < a < 2*N - b/2:
				res.append(HexaPoint(a - 1, b + 1))
			if 0 <= a < 2*N - b/2 - 1:
				res.append(HexaPoint(a, b + 1))
			if 0 <= a < 2*N - b/2:
				res.append(HexaPoint(a, b - 1))
		elif b-1 in np.arange(0, 2*N, 2):
			if 0 <= a < 2*N - b/2 - 1:
				res.append(HexaPoint(a, b - 1))
				res.append(HexaPoint(a + 1, b - 1))
			if (b < 2*N - 1) and (0<= a < 2*N - b/2 - 1):
				res.append(HexaPoint(a, b + 1))
		elif -b in np.arange(2, 2*N + 2, 2):
			if 0 <= a < 2*N + b/2:
				res.append(HexaPoint(a, b + 1))
				res.append(HexaPoint(a + 1, b + 1))
			if (b > -2*N) and (0 <= a < 2*N + b/2):
				res.append(HexaPoint(a, b - 1))
		elif -b-1 in np.arange(0, 2*N, 2):
			if 0 <= a < 2*N + b/2:
				res.append(HexaPoint(a, b + 1))
			if 0 <= a < 2*N + b/2 - 1:
				res.append(HexaPoint(a, b - 1))
			if 0 < a < 2*N + b/2:
				res.append(HexaPoint(a - 1, b - 1))
		return res



class Dimer:
	"""
	A container for the dimers (pairs of neighbors in a honeycomb).

	Useful methods:
	- draw_on_ax(ax)
	- type(): either "up", "left" or "right".
	"""
	def __init__(self, p, q, N):
		self.link = []
		if q in p.neighbors(N):
			if p.y < q.y:
				self.link = [p, q]
			else:
				self.link = [q, p]

	def __eq__(self, other):
		return (self.link == other.link)

	def draw_on_ax(self, ax, beautiful=True):
		L = self.link
		if len(L) == 2:
			if beautiful:
				x = L[0].graphic_coordinates()[0]
				y = L[0].graphic_coordinates()[1]
				if self.type() == "up":
					Lx = [x + 1/(2*np.sqrt(3)), x + 1/(2*np.sqrt(3)), x - 1/(np.sqrt(3)), x - 1/(np.sqrt(3))]
					Ly = [y - 1/2, y + 1/2, y+1, y]
					ax.fill(Lx, Ly, color=(0.8,0.8,0.8))
				elif self.type() == "left":
					Lx = [x + 1/(np.sqrt(3)), x - 1/(2*np.sqrt(3)), x - 2/np.sqrt(3), x - 1/(2*np.sqrt(3))]
					Ly = [y, y + 1/2, y, y - 1/2]
					ax.fill(Lx, Ly, color=(0.7,0.9,0.9))
				elif self.type() == "right":
					Lx = [x - 1/(2*np.sqrt(3)), x - 1/(2*np.sqrt(3)), x + 1/np.sqrt(3), x + 1/np.sqrt(3)]
					Ly = [y - 1/2, y + 1/2, y+1, y]
					ax.fill(Lx, Ly, color=(0.95,0.95,0.95))
				ax.plot(Lx + [Lx[0]], Ly + [Ly[0]], color="k")
			else:
				p = L[0].graphic_coordinates()
				q = L[1].graphic_coordinates()
				ax.plot([p[0], q[0]], [p[1], q[1]], color="k")

	def type(self):
		xp = (self.link[0]).graphic_coordinates(rotate=False)[0]
		xq = (self.link[1]).graphic_coordinates(rotate=False)[0]
		if xp == xq:
			return "up"
		elif xp < xq:
			return "right"
		else:
			return "left"



class Hexagon:
	"""
	A container for the hexagons that constitute a honeycomb.

	Useful methods:
	- draw_on_ax(ax)
	- points(): the list of 6 points of the hexagon.
	- inner_dimers(N): the 3 dimers of the hexagon when it is an inner corner
	of a dimer covering of the honeycomb of size N.
	- outer_dimers(N): the 3 dimers of the hexagon when it is an outer corner
	of a dimer covering of the honeycomb of size N.
	"""
	def __init__(self, a, b):
		self.x = a
		self.y = b

	def __repr__(self):
		return "H" + str([self.x, int(self.y)])

	def __eq__(self, other):
		return ((self.x == other.x) and (self.y == other.y))

	def draw_on_ax(self, ax):
		L = [p.graphic_coordinates() for p in self.points()]
		Lx, Ly = [p[0] for p in L] + [L[0][0]], [p[1] for p in L] + [L[0][1]]
		ax.plot(Lx, Ly, color="k")
		
	def points(self):
		a = self.x
		b = self.y
		if b > 0:
			return [HexaPoint(a, 2*b), HexaPoint(a, 2*b + 1), HexaPoint(a + 1, 2*b),
				HexaPoint(a + 1, 2*b - 1), HexaPoint(a + 1, 2*b - 2), HexaPoint(a, 2*b - 1)]
		elif b == 0:
			return [HexaPoint(a, 2*b), HexaPoint(a, 2*b + 1), HexaPoint(a + 1, 2*b),
				HexaPoint(a + 1, 2*b - 1), HexaPoint(a, 2*b - 2), HexaPoint(a, 2*b - 1)]
		elif b < 0:
			return [HexaPoint(a, 2*b), HexaPoint(a + 1, 2*b + 1), HexaPoint(a + 1, 2*b),
				HexaPoint(a + 1, 2*b - 1), HexaPoint(a, 2*b - 2), HexaPoint(a, 2*b - 1)]

	def inner_dimers(self, N):
		L = self.points()
		res = [Dimer(L[0], L[1], N),
			   Dimer(L[2], L[3], N),
			   Dimer(L[4], L[5], N)]
		return res

	def outer_dimers(self, N):
		L = self.points()
		res = [Dimer(L[1], L[2], N),
			   Dimer(L[3], L[4], N),
			   Dimer(L[5], L[0], N)]
		return res

	def exists(self):
		return not (self.x == -1)

	def numbered_neighbor(self, N, k=1):
		a = self.x
		b = self.y
		if k == 1:
			if ((0 <= b <= N - 1) and (0 < a < 2*N - b - 1)):
				return Hexagon(a - 1, b)
			elif ((-N + 1 <= b < 0) and (0 < a < 2*N + b - 1)):
				return Hexagon(a - 1, b)
			else:
				return Hexagon(-1, 0)
		elif k == 2:
			if ((0 <= b < N - 1) and (0 < a < 2*N - b - 1)):
				return Hexagon(a - 1, b + 1)
			elif ((-N + 1 <= b < 0) and (0 <= a < 2*N + b - 1)):
				return Hexagon(a, b + 1)
			else:
				return Hexagon(-1, 0)
		elif k == 3:
			if ((0 <= b < N - 1) and (0 <= a < 2*N - b - 2)):
				return Hexagon(a, b + 1)
			elif ((-N + 1 <= b < 0) and (0 <= a < 2*N + b - 1)):
				return Hexagon(a + 1, b + 1)
			else:
				return Hexagon(-1, 0)
		elif k == 4:
			if ((0 <= b <= N - 1) and (0 <= a < 2*N - b - 2)):
				return Hexagon(a + 1, b)
			elif ((-N + 1 <= b < 0) and (0 <= a < 2*N + b - 2)):
				return Hexagon(a + 1, b)
			else:
				return Hexagon(-1, 0)
		elif k == 5:
			if ((0 < b <= N - 1) and (0 <= a < 2*N - b - 1)):
				return Hexagon(a + 1, b - 1)
			elif ((-N + 1 < b <= 0) and (0 <= a < 2*N + b - 2)):
				return Hexagon(a, b - 1)
			else:
				return Hexagon(-1, 0)
		elif k == 6:
			if ((0 < b <= N - 1) and (0 <= a < 2*N - b - 1)):
				return Hexagon(a, b - 1)
			elif ((-N + 1 < b <= 0) and (0 < a < 2*N + b - 1)):
				return Hexagon(a - 1, b - 1)
			else:
				return Hexagon(-1, 0)

	def left_neighbors(self, N):
		res = []
		for k in [1,3,5]:
			h = self.numbered_neighbor(N, k)
			if h.exists():
				res.append(h)
		return res

	def right_neighbors(self, N):
		res = []
		for k in [2,4,6]:
			h = self.numbered_neighbor(N, k)
			if h.exists():
				res.append(h)
		return res



class Honeycomb:
	"""
	Honeycomb(N) returns the honeycomb of size N. Useful methods:
	- draw_on_ax(ax)
	- plot()
	- points(): the list of points of the lattice.
	- hexagons(): the list of hexagons in the lattice.
	- random_hexagon(): a random hexagon in the lattice.
	"""
	def __init__(self, N):
		self.size = N

	def __repr__(self):
		return f"Honeycomb lattice of size {self.size}"

	def draw_on_ax(self, ax, with_hexagons=True, with_text=False):
		if with_hexagons:
			for H in self.hexagons():
				H.draw_on_ax(ax)
		else:
			L = [p.graphic_coordinates() for p in self.points()]
			Lx, Ly = [p[0] for p in L], [p[1] for p in L]
			ax.scatter(Lx, Ly, color="k")
		if with_text:
			for H in self.hexagons():
				coor = HexaPoint(H.x, 2*H.y).graphic_coordinates()
				ax.text(coor[0] + 1/(2*np.sqrt(3)), coor[1], str((H.x, H.y)))

	def plot(self, with_hexagons=True, with_text=False):
		fig, ax = plt.subplots(figsize=(10, 10))
		self.draw_on_ax(ax, with_hexagons, with_text)
		ax.set_aspect(1)
		ax.set_axis_off()
		plt.show()

	def points(self):
		res = []
		N = self.size
		for j in range(N):
			for i in range(0, 2*N - j):
				res.append(HexaPoint(i, 2*j))
				res.append(HexaPoint(i, -2*j - 1))
			for i in range(0, 2*N - j - 1):
				res.append(HexaPoint(i, 2*j + 1))
				res.append(HexaPoint(i, -2*j - 2))
		return res

	def hexagons(self):
		res = []
		N = self.size
		for j in range(N):
			for i in range(2*N - j - 1):
				res.append(Hexagon(i, j))
		for j in range(1, N):
			for i in range(2*N - j - 1):
				res.append(Hexagon(i, -j))
		return res

	def random_hexagon(self):
		L = self.hexagons()
		nb = len(L)
		return L[rand.randint(nb)]





class DimerCovering:
	"""
	A container for dimer coverings of a honeycomb lattice. DimerCovering(N)
	returns the empty dimer covering of the honeycomb of size N, whereas
	DimerCovering(N, full=True) returns the full dimer covering.

	Useful methods:
	- plot()
	- plot_honeycomb()
	- plane_partition()
	- volume()
	- switchup(H): add the box corresponding to the hexagon H, if possible.
	- switchdown(H): remove the box corresponding to the hexagon H, if possible.
	- random_local_switch(H, q): performs a random switch at H, driven by a Bernoulli
	 random variable of parameter q/(q+1).
	- transform(n, q): performs n transitions of the Markov chain with parameter q on
	the set of all dimer coverings with same size.
	"""
	def __init__(self, N, full=False):
		self.size = N
		self.dimers = []
		if full:
			self.inner_corners = {}
			self.outer_corners = {}
			self.grid = N * np.ones((N, N))
			self.outer_corners[(N - 1, 0)] = (N - 1, N - 1)
			for j in range(N):
				for i in range(N - j - 1,2*N - j - 1):
					p = HexaPoint(i, 2*j + 1)
					q = HexaPoint(i + 1, 2*j)
					self.dimers.append(Dimer(p, q, N))
					p = HexaPoint(i, -2*j - 2)
					q = HexaPoint(i + 1, -2*j - 1)
					self.dimers.append(Dimer(p, q, N))
				for i in range(N - j):
					p = HexaPoint(i, 2*j)
					q = HexaPoint(i, 2*j - 1)
					self.dimers.append(Dimer(p, q, N))
			for j in range(1, N):
				for i in range(N - j):
					p = HexaPoint(i, -2*j)
					q = HexaPoint(i, -2*j - 1)
					self.dimers.append(Dimer(p, q, N))
		else:
			self.inner_corners = {}
			self.outer_corners = {}
			self.grid = np.zeros((N, N))
			self.inner_corners[(N - 1, 0)] = (0,0)
			for i in range(N):
				for j in range(N):
					p = HexaPoint(i, 2*j)
					q = HexaPoint(i, 2*j + 1)
					self.dimers.append(Dimer(p, q, N))
					p = HexaPoint(i, -2*j - 1)
					q = HexaPoint(i, -2*j - 2)
					self.dimers.append(Dimer(p, q, N))
			for j in range(N):
				for i in range(N, 2*N - j):
					p = HexaPoint(i, 2*j)
					q = HexaPoint(i, 2*j - 1)
					self.dimers.append(Dimer(p, q, N))
			for j in range(1, N):
				for i in range(N, 2*N - j):
					p = HexaPoint(i, -2*j)
					q = HexaPoint(i, -2*j - 1)
					self.dimers.append(Dimer(p, q, N))

	def __repr__(self):
		return f"Dimer covering of the honeycomb of size {self.size}"

	def __eq__(self, other):
		return np.all(self.grid == other.grid)

	def plot(self, beautiful=True):
		fig, ax = plt.subplots(figsize=(10, 10))
		if (not beautiful):
			Honeycomb(self.size).draw_on_ax(ax, with_hexagons=False)
		for D in self.dimers:
			D.draw_on_ax(ax, beautiful)
		ax.set_aspect(1)
		ax.set_axis_off()
		plt.show()		
		
	def plot_honeycomb(self):
		Honeycomb(self.size).plot(with_text=True)

	def volume(self):
		return int(np.sum(self.grid))

	def plane_partition(self):
		return self.grid
		
	def has_for_inner_corner(self, H):
		return all([(D in self.dimers) for D in H.inner_dimers(self.size)])

	def has_for_outer_corner(self, H):
		return all([(D in self.dimers) for D in H.outer_dimers(self.size)])

	def switchup(self, H):
		Hc = (H.x, H.y)
		if Hc in self.inner_corners.keys():
			for J in H.right_neighbors(self.size):
				if self.has_for_outer_corner(J):
					del self.outer_corners[(J.x, J.y)]
			for D in H.inner_dimers(self.size):
				self.dimers.remove(D)
			for D in H.outer_dimers(self.size):
				self.dimers.append(D)
			v = self.inner_corners[Hc]
			i = v[0]
			j = v[1]
			self.grid[i,j] +=1
			del self.inner_corners[Hc]
			self.outer_corners[Hc] = (i, j)
			J1 = H.numbered_neighbor(self.size, 1)
			if J1.exists() and self.has_for_inner_corner(J1):
				self.inner_corners[(J1.x, J1.y)] = (i+1, j)
			J3 = H.numbered_neighbor(self.size, 3)
			if J3.exists() and self.has_for_inner_corner(J3):
				self.inner_corners[(J3.x, J3.y)] = (i, j)
			J5 = H.numbered_neighbor(self.size, 5)
			if J5.exists() and self.has_for_inner_corner(J5):
				self.inner_corners[(J5.x, J5.y)] = (i, j+1)
		else:
			pass

	def switchdown(self, H):
		Hc = (H.x, H.y)
		if Hc in self.outer_corners.keys():
			for J in H.left_neighbors(self.size):
				if self.has_for_inner_corner(J):
					del self.inner_corners[(J.x, J.y)]
			for D in H.outer_dimers(self.size):
				self.dimers.remove(D)
			for D in H.inner_dimers(self.size):
				self.dimers.append(D)
			v = self.outer_corners[Hc]
			i = v[0]
			j = v[1]
			self.grid[i,j] -= 1
			del self.outer_corners[Hc]
			self.inner_corners[Hc] = (i, j)
			J2 = H.numbered_neighbor(self.size, 2)
			if J2.exists() and self.has_for_outer_corner(J2):
				self.outer_corners[(J2.x, J2.y)] = (i, j-1)
			J4 = H.numbered_neighbor(self.size, 4)
			if J4.exists() and self.has_for_outer_corner(J4):
				self.outer_corners[(J4.x, J4.y)] = (i-1, j)
			J6 = H.numbered_neighbor(self.size, 6)
			if J6.exists() and self.has_for_outer_corner(J6):
				self.outer_corners[(J6.x, J6.y)] = (i, j)
		else:
			pass
		
	def random_local_switch(self, H, q=1):
		alea = bool(scs.bernoulli(q/(q+1)).rvs())
		if alea:
			self.switchup(H)
		else:
			self.switchdown(H)

	def transform(self, n=1, q=1):
		Hexagons = Honeycomb(self.size).hexagons()
		for k in range(n):
			H = rand.choice(Hexagons)
			self.random_local_switch(H, q)

	def transform_from_list(self, L, M):
		for k in range(len(L)-1, -1, -1):
			H = L[k]
			if M[k]:
				self.switchup(H)
			else:
				self.switchdown(H)



def RandomDimerCovering(N, q=1):
	"""
	Computes a random surface inside the cube of edge N, of probability proportional to
	q^volume.

	Parameters:
	- N: the size of the cube containing the random surface.
	- q: the parameter of the probability measure.
	"""
	count = 0
	LL = []
	M = []
	Hexagons = Honeycomb(N).hexagons()
	Hlen = len(Hexagons)
	EP = DimerCovering(N, full=False)
	FP = DimerCovering(N, full=True)
	test = True
	while test:
		LL += [rand.choice(Hlen) for k in range(2**count)]
		M += [bool(scs.bernoulli(q/(q+1)).rvs()) for k in range(2**count)]
		count += 1
		EP.transform_from_list([Hexagons[i] for i in LL], M)
		FP.transform_from_list([Hexagons[i] for i in LL], M)
		if EP == FP:
			test = False
		else:
			EP = DimerCovering(N, full=False)
			FP = DimerCovering(N, full=True)
	return EP