Module ateams.models
Models
A model programmatically specifies a probability distribution over a combinatorial object like a graph, cubical complex, or polytope. Given one of these objects (e.g. Cubical), each model implements a proposal function for sampling from the distribution it specifies. For example, SwendsenWang represents the Edwards-Sokal coupling of the plaquette random cluster model and Potts lattice gauge theory on subcomplexes of the $d$-fold, scale-$N$ cubical torus $\mathbb T^d_N$; its proposal function is a generalized version of the Swendsen-Wang algorithm.
The SwendsenWang, InvadedCluster, and Nienhuis models require performant sparse linear algebra over finite fields, for which we turn to LinBox. The InvadedCluster and Bernoulli models require fast computation of persistent homology, for which we use PHAT (when sampling over $\mathbb Z/2\mathbb Z$) and a home-brew variant of the "twist" algorithm proposed by Chen and Kerber. Only the Glauber dynamics model is implemented in "pure Python," though it relies exclusively on NumPy.
The documentation for each model comes with performance information, viewable by clicking (e.g.) the "Performance over $\mathbb T^2_N$" button. Each configuration completed 100 iterations on Pangolin, a Dell Precision 5280 workstation with an 18-core Intel Xeon W-2295 CPU clocked at 1.3GHz.
Classes
class Bernoulli (C, dimension=1, **kwargs)-
Initializes classic Bernoulli percolation on the provided complex, detecting percolation in the
dimension-1th homology group.Args
C:Complex- The
Complexobject on which we'll be running experiments. dimension:int=1- The dimension of cells on which we're percolating.
Expand source code
class Bernoulli(): _name = "Bernoulli" def __init__(self, C, dimension=1, **kwargs): """ Initializes classic Bernoulli percolation on the provided complex, detecting percolation in the `dimension`-1th homology group. Args: C (Complex): The `Complex` object on which we'll be running experiments. dimension (int=1): The dimension of cells on which we're percolating. """ # Object access. Set a field. self.complex = C self.dimension = dimension self._returns = 1 self.field = 2 # Set a phantom spins attribute so we don't break the Chain self.spins = None # Force-recompute the matrices for a different dimension; creates # a set of orientations for fast elementwise products. self.matrices = Matrices() self.matrices.full = self.complex.matrices.full boundary, coboundary = self.complex.recomputeBoundaryMatrices(dimension) self.matrices.boundary = boundary self.matrices.coboundary = coboundary # Useful values to have later. self.cellCount = len(self.complex.flattened) self.cells = len(self.complex.Boundary[self.dimension]) self.faces = len(self.complex.Boundary[self.dimension-1]) self.target = np.arange(self.complex.breaks[self.dimension], self.complex.breaks[self.dimension+1]) # Premake the "occupied cells" array; change the dimension of the complex # to correspond to the provided dimension. self.rank = comb(len(self.complex.corners), self.dimension) self.nullity = len(self.complex.Boundary[self.dimension]) # Delegates computation for persistence. self._delegateComputation() # Seed the random number generator. self.RNG = np.random.default_rng() def _delegateComputation(self): low, high = self.complex.breaks[self.dimension], self.complex.breaks[self.dimension+1] times = set(range(self.cellCount)) def whittle(pairs): _births, _deaths = zip(*pairs) births = set(_births) deaths = set(_deaths) return set( e for e in times-(births|deaths) if low <= e < high ) def persist(filtration): essential = ComputePersistencePairs( self.matrices.full, filtration, self.dimension, self.complex.breaks ) return whittle(essential) self.persist = persist def _filtrate(self): """ Constructs a filtration based on the evaluation of the cochain. """ # Find which cubes get zeroed out (i.e. are sent to zero by the cocycle). low = self.complex.breaks[self.dimension] high = self.complex.breaks[self.dimension+1] filtration = np.arange(self.cellCount) shuffled = np.random.permutation(self.target) filtration[low:high] = shuffled return filtration, shuffled-low def _initial(self): """ Computes an initial state for the model's Complex. Returns: A numpy `np.array` representing a vector of spin assignments. """ return self.RNG.integers( 0, high=self.field, dtype=FINT, size=self.faces ) def _proposal(self, time): """ Proposal scheme for generalized invaded-cluster evolution on the random-cluster model. Args: time (int): Step in the chain. Returns: A numpy array representing a vector of spin assignments. """ # Construct the filtration and find the essential cycles. filtration, shuffled = self._filtrate() essential = self.persist(filtration) j = 0 low = self.complex.breaks[self.dimension] occupied = np.zeros((self.rank, self.nullity)) for t in essential: occupiedIndices = shuffled[:t-low] occupied[j,occupiedIndices] = 1 j += 1 return occupied def _assign(self, cocycle): pass class Glauber (C,
dimension=1,
field=2,
temperature=<function constant.<locals>._>,
initial=None,
**kwargs)-
Initializes Glauber dynamics on the Potts model.
Args
C- The
Complexobject on which we'll be running experiments. dimension:int=1- Dimension on which we'll be building experiments.
field:int=2- Field characteristic.
temperature:Callable- A temperature schedule function which
takes a single positive integer argument
t, and returns the scheduled temperature at timet. initial:np.ndarray- A vector of spin assignments to components.
s/iteration $N$ Cells Field 256 131,072 $\mathbb Z/2\mathbb Z$ 0.007 $\mathbb Z/3\mathbb Z$ 0.008 $\mathbb Z/5\mathbb Z$ 0.007 $\mathbb Z/7\mathbb Z$ 0.007 512 524,288 $\mathbb Z/2\mathbb Z$ 0.056 $\mathbb Z/3\mathbb Z$ 0.023 $\mathbb Z/5\mathbb Z$ 0.026 $\mathbb Z/7\mathbb Z$ 0.029 1024 2,097,152 $\mathbb Z/2\mathbb Z$ 0.234 $\mathbb Z/3\mathbb Z$ 0.108 $\mathbb Z/5\mathbb Z$ 0.104 $\mathbb Z/7\mathbb Z$ 0.110 s/iteration $N$ Cells Field 20 960,000 $\mathbb Z/2\mathbb Z$ 0.078 $\mathbb Z/3\mathbb Z$ 0.065 $\mathbb Z/5\mathbb Z$ 0.071 $\mathbb Z/7\mathbb Z$ 0.078 30 4,860,000 $\mathbb Z/2\mathbb Z$ 0.376 $\mathbb Z/3\mathbb Z$ 0.379 $\mathbb Z/5\mathbb Z$ 0.383 $\mathbb Z/7\mathbb Z$ 0.367 40 15,360,000 $\mathbb Z/2\mathbb Z$ 1.271 $\mathbb Z/3\mathbb Z$ 1.172 $\mathbb Z/5\mathbb Z$ 1.147 $\mathbb Z/7\mathbb Z$ 1.185 Expand source code
class Glauber(): _name = "Glauber" def __init__( self, C, dimension=1, field=2, temperature=constant(-0.6), initial=None, **kwargs ): """ Initializes Glauber dynamics on the Potts model. Args: C: The `Complex` object on which we'll be running experiments. dimension (int=1): Dimension on which we'll be building experiments. field (int=2): Field characteristic. temperature (Callable): A temperature schedule function which takes a single positive integer argument `t`, and returns the scheduled temperature at time `t`. initial (np.ndarray): A vector of spin assignments to components. <center> <button type="button" class="collapsible" id="Glauber-_proposal-2"> Performance in \(\mathbb T^2_N\)</button> </center> ..include:: ./tables/Glauber._proposal.2.html <center> <button type="button" class="collapsible" id="Glauber-_proposal-4"> Performance in \(\mathbb T^4_N\)</button> </center> ..include:: ./tables/Glauber._proposal.4.html """ self.field = field self.complex = C self.dimension = dimension self.temperature = temperature self._returns = 2 # Useful values to have later. self.cells = len(self.complex.Boundary[self.dimension]) self.faces = len(self.complex.Boundary[self.dimension-1]) self.closed = 1 # Seed the random number generator. self.RNG = np.random.default_rng() # If no initial spin configuration is passed, initialize. if not initial: self.spins = self._initial() else: self.spins = (initial%self.field).astype(FINT) def _initial(self): """ Computes an initial state for the model's Complex. Returns: A numpy `np.array` representing a vector of spin assignments. """ return self.RNG.integers( 0, high=self.field, dtype=MINT, size=self.faces ) def _proposal(self, time): """ Proposal scheme for generalized Glauber dynamics on the Potts model: uniformly randomly chooses a face in the complex, flips the face's spin, and returns the corresponding cocycle. Args: time (int): Step in the chain. Returns: A NumPy array representing a vector of spin assignments. """ # Flip. spins = self.spins[::] flip = self.RNG.integers(self.faces) spins[flip] = (self.RNG.integers(self.field, dtype=FINT)) boundary = self.complex.Boundary[self.dimension] q = self.field # Evaluate the current spin assignment (cochain). evaluation = spins[boundary] evaluation[:, 1::2] = -evaluation[:, 1::2]%q sums = evaluation.sum(axis=1)%q opens = np.nonzero(sums == 0)[0] closed = -(self.faces-opens.shape[0]) # If we don't know how many "closed" faces there are, compute it. if self.closed > 0: evaluation = self.spins[boundary] evaluation[:,1::2] = -evaluation[:,1::2]%q sums = evaluation.sum(axis=1)%q opens = np.nonzero(sums == 0)[0] self.closed = -(self.faces-opens.shape[0]) # Compute the energy and decide if this is an acceptable transition. energy = np.exp(self.temperature(time)*(self.closed-closed)) if self.RNG.uniform() < min(1, energy): self.closed = closed self.open = opens self.spins = spins satisfied = np.zeros(self.cells, dtype=FINT) satisfied[self.open] = 1 return self.spins, satisfied def _assign(self, cochain): """ Updates mappings from faces to spins and cubes to occupations. Args: cochain (np.array): Cochain on the complex. Returns: None. """ self.spins = cochain class InvadedCluster (C,
dimension=1,
field=2,
initial=None,
full=False,
stop=<function InvadedCluster.<lambda>>,
tries=64,
**kwargs)-
Initializes the plaquette invaded-cluster algorithm on the provided integer complex, detecting percolation in the
homology-th homology group.Args
C:Complex- The
Complexobject on which we'll be running experiments. dimension:int=1- The dimension of cells on which we're percolating.
field:int=2- Field characteristic.
initial:np.array- A vector of spin assignments to components.
full:bool=True- Do we find all the homological persistence events, or just the target one?
stop:function- A function that returns the number of essential cycles found before sampling the next configuration.
tries:int=64- Number of persistence computation attempts made before re-sampling.
s/iteration $N$ Cells Field 256 131,072 $\mathbb Z/2\mathbb Z$ 0.237 $\mathbb Z/3\mathbb Z$ 12.595 $\mathbb Z/5\mathbb Z$ 8.843 $\mathbb Z/7\mathbb Z$ 7.768 512 524,288 $\mathbb Z/2\mathbb Z$ 1.105 $\mathbb Z/3\mathbb Z$ 145.532 $\mathbb Z/5\mathbb Z$ 76.667 $\mathbb Z/7\mathbb Z$ 62.254 s/iteration $N$ Cells Field 10 60,000 $\mathbb Z/2\mathbb Z$ 0.847 $\mathbb Z/3\mathbb Z$ 42.541 Expand source code
class InvadedCluster(): _name = "InvadedCluster" def __init__( self, C, dimension=1, field=2, initial=None, full=False, stop=lambda: 1, tries=64, **kwargs ): """ Initializes the plaquette invaded-cluster algorithm on the provided integer complex, detecting percolation in the `homology`-th homology group. Args: C (Complex): The `Complex` object on which we'll be running experiments. dimension (int=1): The dimension of cells on which we're percolating. field (int=2): Field characteristic. initial (np.array): A vector of spin assignments to components. full (bool=True): Do we find *all* the homological persistence events, or just the target one? stop (function): A function that returns the number of essential cycles found before sampling the next configuration. tries (int=64): Number of persistence computation attempts made before re-sampling. <center> <button type="button" class="collapsible" id="InvadedCluster-Persistence-2">Performance in \(\mathbb T^2_N\)</button> </center> ..include:: ./tables/InvadedCluster.Persistence.2.html <center> <button type="button" class="collapsible" id="InvadedCluster-Persistence-4">Performance in \(\mathbb T^4_N\)</button> </center> ..include:: ./tables/InvadedCluster.Persistence.4.html """ # Object access. self.complex = C self.dimension = dimension self.stop = stop self._returns = 3 self.field = field self.tries = tries self.full = full # Check if we're debugging. self._DEBUG = kwargs.get("_DEBUG", False) # Force-recompute the matrices for a different dimension; creates # a set of orientations for fast elementwise products. self.matrices = Matrices() self.matrices.full = self.complex.matrices.full boundary, coboundary = self.complex.recomputeBoundaryMatrices(dimension) self.matrices.boundary = boundary self.matrices.coboundary = coboundary # Useful values to have later. self.cellCount = len(self.complex.flattened) self.cells = len(self.complex.Boundary[self.dimension]) self.faces = len(self.complex.Boundary[self.dimension-1]) self.target = np.arange(self.complex.breaks[self.dimension], self.complex.breaks[self.dimension+1]) # Check the dimensions of the boundary/coboundary matrices by comparing # the number of cells. LinBox is really sensitive to smaller-size matrices, # but can easily handle large ones. if self.cells*self.faces < 100000: warnings.warn(f"complex with {self.cells*self.faces} boundary matrix entries is too small for accurate matrix solves; may segfault.", TooSmallWarning, stacklevel=2) # Premake the "occupied cells" array; change the dimension of the complex # to correspond to the provided dimension. self.rank = comb(len(self.complex.corners), self.dimension) # Delegates computation for persistence and cocycle sampling. self._delegateComputation() # Seed the random number generator. self.RNG = np.random.default_rng() # If no initial spin configuration is passed, initialize. if not initial: self.spins = self._initial() else: self.spins = (initial%self.field).astype(FINT) def _delegateComputation(self): low, high = self.complex.breaks[self.dimension], self.complex.breaks[self.dimension+1] # If we're using LinBox, our sampling method is Lanczos regardless of # dimension. def sample(zeros): try: return np.array(ReducedKernelSample( self.matrices.coboundary, zeros, 2*self.dimension, self.faces, self.field, self._DEBUG ), dtype=FINT) except Exception as e: raise NumericalInstabilityWarning(e) # If we're using LinBox and the characteristic of our field is greater # than 2, we use the twist_reduce variant implemented in this library. if self.field > 2: Twister = Twist(self.field, self.matrices.full, self.complex.breaks, self.cellCount, self.dimension, self._DEBUG) Twister.LinearComputeCobasis(); def _user_continue_prompt(): _affirm = {"", "yes", "y", "ye"} _cont = input("\ncontinue? [y/n, default y] --> ") if _cont not in _affirm: sys.exit(1) print("\n\n") def persist(filtration): # Set defaults. essential = set() tries = 0 # Debugging mode. Compares the SparseRREF/SpaSM output to the # classical persistence algorithm (LinearComputePercolationEvents). # NOTE: this waits for user input to continue after each iteration. if self._DEBUG: print("######################################") print("#### SPARSERREF/SPASM COMPUTATION ####") print("######################################") rstart = time.time() att = 0 while len(essential) < self.rank: stop = 0 if self.full else self._STOP essential = Twister.RankComputePercolationEvents(filtration, stop=stop) att += 1 if not self.full and len(essential) == 1: break; rend = time.time() essential = np.array(list(essential)) essential = essential[(essential >= low) & (essential < high)] essential.sort() print("\n\n") print(essential) print(f"time: {(rend-rstart):.4f} :: {att} attempts") _user_continue_prompt() # print("######################################") # print("#### LINEAR COMPUTATION ####") # print("######################################") # lstart = time.time() # _essential = Twister.LinearComputePercolationEvents(filtration) # lend = time.time() # _essential = np.array(list(_essential)) # _essential = _essential[(_essential >= low) & (_essential < high)] # _essential.sort() # print("\n\n") # print(_essential, f"time: {(lend-lstart):.4f}") # _user_continue_prompt() else: # Attempt the persistence computation up to `self.tries` times; # if we exceed the attempts budget, re-sample the filtration and # try again. while len(essential) < self.rank: # Compute essential cycle birth times. If we aren't computing the # stop = 0 if self.full else self._STOP essential = Twister.RankComputePercolationEvents(filtration, stop=stop) tries += 1 # If we don't have enough essential cycles and we've tried # too many times, re-sample the filtration and try again. if len(essential) < self.rank and tries >= self.tries: print(f"[Persistence] exceeded the acceptable number of {self.tries} attempts. Re-sampling filtration...") filtration = self._filtrate(self.spins) tries = 0 if not self.full and len(essential) == 1: break essential = np.array(list(essential)) essential = essential[(essential >= low) & (essential < high)] essential.sort() return essential # If we're using LinBox and the field we're computing over *is* two, # use PHAT. elif self.field < 3: times = set(range(self.cellCount)) def whittle(pairs): _births, _deaths = zip(*pairs) births = set(_births) deaths = set(_deaths) return set( e for e in times-(births|deaths) if low <= e < high ) def persist(filtration): essential = ComputePersistencePairs(self.matrices.full, filtration, self.dimension, self.complex.breaks) return whittle(essential) self.sample = sample self.persist = persist def _filtrate(self, cochain): """ Constructs a filtration based on the evaluation of the cochain. """ # Find which cubes get zeroed out (i.e. are sent to zero by the cocycle). boundary = self.complex.Boundary[self.dimension] q = self.field coefficients = cochain[boundary] coefficients[:,1::2] = -coefficients[:,1::2]%q sums = coefficients.sum(axis=1)%q # POSSIBLY INEFFICIENT!! TAKE A LOOK DUMMY satisfied = np.nonzero(sums==0)[0] unsatisfied = np.nonzero(sums>0)[0] m = satisfied.shape[0] # Construct the filtration. filtration = np.arange(self.cellCount) low = self.complex.breaks[self.dimension] high = self.complex.breaks[self.dimension+1] shuffled = np.random.permutation(satisfied) filtration[low:low+m] = self.target[shuffled] filtration[low+m:high] = self.target[unsatisfied] return filtration, shuffled, satisfied def _initial(self): """ Computes an initial state for the model's Complex. Returns: A numpy `np.array` representing a vector of spin assignments. """ return self.RNG.integers( 0, high=self.field, dtype=FINT, size=self.faces ) def _proposal(self, time): """ Proposal scheme for generalized invaded-cluster evolution on the random-cluster model. Args: time (int): Step in the chain. Returns: A numpy array representing a vector of spin assignments. """ # Set a stopping time. self._STOP = self.stop(); # Construct the filtration and find the essential cycles. filtration, shuffledIndices, satisfiedIndices = self._filtrate(self.spins) essential = self.persist(filtration) j = 0 low = self.complex.breaks[self.dimension] satisfied = np.zeros(self.cells, dtype=int) # If we only computed a desired percolation time, then the matrix of # occupied cell indices has one row... if not self.full: t = essential[0] occupied = np.zeros(self.cells, dtype=int) occupiedIndices = shuffledIndices[:t-low+1] occupied[occupiedIndices] = 1 spins = self.sample(occupiedIndices) # Otherwise, it has `rank` rows. else: occupied = np.zeros((self.rank, self.cells), dtype=int) for t in sorted(essential): occupiedIndices = shuffledIndices[:t-low+1] occupied[j,occupiedIndices] = 1 if (j+1) == self._STOP: spins = self.sample(occupiedIndices) j += 1 satisfied[satisfiedIndices] = 1 return spins, occupied, satisfied def _assign(self, cocycle): """ Updates mappings from faces to spins and cubes to occupations. Args: cocycle (np.array): Cocycle on the subcomplex. Returns: Nothing. """ self.spins = cocycle class Nienhuis (C, q1, q2, dimension=2, field=2, initial=None, maxTries=16, **kwargs)-
Initializes the self-dual Nienhuis model.
Args
C:Complex- The
Complexobject on which we'll be running experiments. Here, the dimension is assumed to be 2. q1:float- Edge coupling parameter.
q2:float- Face (i.e. cell) coupling parameter.
dimension:int- Dimension targeted by the model.
field:int=2- Field characteristic.
initial:np.array- A vector of spin assignments to components.
maxTries:int=16- The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix.
Expand source code
class Nienhuis(): _name = "Nienhuis" def __init__( self, C, q1, q2, dimension=2, field=2, initial=None, maxTries=16, **kwargs ): """ Initializes the self-dual Nienhuis model. Args: C (Complex): The `Complex` object on which we'll be running experiments. Here, the dimension is assumed to be 2. q1 (float): Edge coupling parameter. q2 (float): Face (i.e. cell) coupling parameter. dimension (int): Dimension targeted by the model. field (int=2): Field characteristic. initial (np.array): A vector of spin assignments to components. maxTries (int=16): The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix. """ # Object access. self.complex = C self.dimension = dimension self._returns = 3 self.field = field # Force-recompute the matrices for a different dimension; creates # a set of orientations for fast elementwise products. Here, the dimension # is implicitly assumed to be 2. self.matrices = Matrices() self.matrices.full = self.complex.matrices.full boundary, coboundary = self.complex.recomputeBoundaryMatrices(self.dimension) self.matrices.boundary = boundary self.matrices.coboundary = coboundary # Useful values to have later. self.cells = len(self.complex.Boundary[self.dimension]) self.faces = len(self.complex.Boundary[self.dimension-1]) self.faceIndices = np.arange(self.faces) self.orientations = Bunch() self.orientations.cells = np.tile([-1,1], self.dimension).astype(FINT) self.orientations.faces = np.tile([-1,1], self.dimension-1).astype(FINT) # Coupling parameters. self.coupling = Bunch() self.coupling.face = q1 self.coupling.cell = q2 self.prob = Bunch() self.prob.face = self.coupling.face/(1+self.coupling.face) self.prob.cell = self.coupling.cell/(1+self.coupling.cell) # Check the dimensions of the boundary/coboundary matrices by comparing # the number of cells. LinBox is really sensitive to smaller-size matrices, # but can easily handle large ones. if self.cells*self.faces < 100000: warnings.warn(f"complex with {self.cells*self.faces} boundary matrix entries is too small for accurate matrix solves; may segfault.", TooSmallWarning, stacklevel=2) # Seed the random number generator. self.RNG = np.random.default_rng() # If no initial spin configuration is passed, initialize. if not initial: self.spins = self._initial() else: self.spins = (initial%self.field).astype(FINT) # Delegate computation. self._delegateComputation(maxTries) def _delegateComputation(self, maxTries): def sample(faceZeros, cellZeros): try: return np.array(SubReducedKernelSample( self.matrices.coboundary, cellZeros, faceZeros, self.field, False ), dtype=FINT) except Exception as e: raise NumericalInstabilityWarning(e) self.sample = sample def _initial(self): """ Computes an initial state for the model's Complex. Returns: A numpy array of spin assignments (mod p). """ return self.RNG.integers( 0, high=self.field, dtype=FINT, size=self.faces ) def _proposal(self, time): """ Proposal scheme for the Nienhuis model. Args: time (int): Step in the chain. Returns: A numpy array of spin assignments (mod p). """ # Choose which cells and faces to include. cells = self.RNG.uniform(size=self.cells) faces = self.RNG.uniform(size=self.faces) threshCells = (cells < self.prob.cell).nonzero()[0] threshFaces = (faces < self.prob.face).nonzero()[0] q = self.field # Evaluate the current spin configuration (cochain) on the boundary # of our chosen cells, and find which to include. cellBoundary = self.complex.Boundary[self.dimension] cellCoefficients = self.spins[cellBoundary] cellCoefficients[:,1::2] = -cellCoefficients[:,1::2]%q cellSums = cellCoefficients.sum(axis=1)%q cellZeros = np.nonzero(cellSums[threshCells] == 0)[0] # Decide which faces (by default, edges) to include, take a submatrix # of the coboundary matrix by including only the rows corresponding to # included cells and columns corresponding to included faces, then # sample new spins. faceCoefficients = np.intersect1d((self.spins==0).nonzero()[0], threshFaces) faceZeros = np.setdiff1d(self.faceIndices, threshFaces) subspins = self.sample(faceZeros, cellZeros) # Organize spin data. spins = self.spins.copy() spins[faceZeros] = subspins spins[faceCoefficients] = 0 # Compute energies. faceEnergy = (self.faces - np.count_nonzero(spins==0))/self.faces cellEnergy = (self.cells - np.count_nonzero(cellSums==0))/self.cells return spins, faceEnergy, cellEnergy def _assign(self, cocycle): """ Updates mappings from faces to spins and cubes to occupations. Args: cocycle (np.array): Cocycle on the subcomplex. Returns: None. """ self.spins = cocycle class SwendsenWang (C,
dimension=1,
field=2,
temperature=<function constant.<locals>._>,
initial=None,
dtypes=(<class 'numpy.int8'>, <class 'numpy.int8'>, <class 'numpy.int8'>),
**kwargs)-
Initializes Swendsen-Wang evolution on the Potts model.
Args
C:Complex- The
Complexobject on which we'll be running experiments. dimension:int=1- The dimension of cells on which we're percolating.
field:int=2- Field characteristic.
temperature:Callable- A temperature schedule function which
takes a single positive integer argument
t, and returns the scheduled temperature at timet. initial:np.array- A vector of spin assignments to components.
dtypes (tuple=(np.int,np.byte,np.byte)): A tuple of NumPy data types to which output data are cast.
s/iteration $N$ Cells Field 256 131,072 $\mathbb Z/2\mathbb Z$ 0.253 $\mathbb Z/3\mathbb Z$ 0.177 $\mathbb Z/5\mathbb Z$ 0.119 $\mathbb Z/7\mathbb Z$ 0.102 512 524,288 $\mathbb Z/2\mathbb Z$ 1.286 $\mathbb Z/3\mathbb Z$ 0.702 $\mathbb Z/5\mathbb Z$ 0.431 $\mathbb Z/7\mathbb Z$ 0.353 1024 2,097,152 $\mathbb Z/2\mathbb Z$ 7.135 $\mathbb Z/3\mathbb Z$ 3.698 $\mathbb Z/5\mathbb Z$ 2.038 $\mathbb Z/7\mathbb Z$ 1.605 2048 8,388,608 $\mathbb Z/2\mathbb Z$ 67.837 $\mathbb Z/3\mathbb Z$ 28.235 $\mathbb Z/5\mathbb Z$ 11.285 $\mathbb Z/7\mathbb Z$ 7.813 s/iteration $N$ Cells Field 10 60,000 $\mathbb Z/2\mathbb Z$ 0.206 $\mathbb Z/3\mathbb Z$ 0.132 $\mathbb Z/5\mathbb Z$ 0.093 $\mathbb Z/7\mathbb Z$ 0.091 20 960,000 $\mathbb Z/2\mathbb Z$ 3.161 $\mathbb Z/3\mathbb Z$ 1.683 $\mathbb Z/5\mathbb Z$ 1.076 $\mathbb Z/7\mathbb Z$ 0.888 30 4,860,000 $\mathbb Z/2\mathbb Z$ 28.550 $\mathbb Z/3\mathbb Z$ 12.539 $\mathbb Z/5\mathbb Z$ 6.381 $\mathbb Z/7\mathbb Z$ 5.234 40 15,360,000 $\mathbb Z/2\mathbb Z$ 224.594 $\mathbb Z/3\mathbb Z$ 79.565 $\mathbb Z/5\mathbb Z$ 27.017 $\mathbb Z/7\mathbb Z$ 19.159 Expand source code
class SwendsenWang(): _name = "SwendsenWang" def __init__( self, C, dimension=1, field=2, temperature=constant(-0.6), initial=None, dtypes=(np.int8,np.byte,np.byte), **kwargs ): r""" Initializes Swendsen-Wang evolution on the Potts model. Args: C (Complex): The `Complex` object on which we'll be running experiments. dimension (int=1): The dimension of cells on which we're percolating. field (int=2): Field characteristic. temperature (Callable): A temperature schedule function which takes a single positive integer argument `t`, and returns the scheduled temperature at time `t`. initial (np.array): A vector of spin assignments to components. dtypes (tuple=(np.int,np.byte,np.byte)): A tuple of NumPy data types to which output data are cast. <center> <button type="button" class="collapsible" id="SwendsenWang-ReducedKernelSample-2"> Performance in \(\mathbb T^2_N\)</button> </center> .. include:: ./tables/SwendsenWang.ReducedKernelSample.2.html <center> <button type="button" class="collapsible" id="SwendsenWang-ReducedKernelSample-4"> Performance in \(\mathbb T^4_N\)</button> </center> .. include:: ./tables/SwendsenWang.ReducedKernelSample.4.html """ # Object access. self.complex = C self.temperature = temperature self.dimension = dimension self._returns = 3 self.field = field self.dtypes = dtypes # Force-recompute the matrices for a different dimension; creates # a set of orientations for fast elementwise products. self.matrices = Matrices() self.matrices.full = self.complex.matrices.full boundary, coboundary = self.complex.recomputeBoundaryMatrices(dimension) self.matrices.boundary = boundary self.matrices.coboundary = coboundary # Useful values to have later. self.cells = len(self.complex.Boundary[self.dimension]) self.faces = len(self.complex.Boundary[self.dimension-1]) # Check the dimensions of the boundary/coboundary matrices by comparing # the number of cells. LinBox is really sensitive to smaller-size matrices, # but can easily handle large ones. if self.cells*self.faces < 10000: warnings.warn(f"complex with {self.cells*self.faces} boundary matrix entries is too small for accurate matrix solves; may segfault.", TooSmallWarning, stacklevel=2) # Seed the random number generator. self.RNG = np.random.default_rng() # Check if we're debugging. self._DEBUG = kwargs.get("_DEBUG", False) # If no initial spin configuration is passed, initialize. if not initial: self.spins = self._initial() else: self.spins = (initial%self.field).astype(FINT) # Delegate computation. self._delegateComputation() def _delegateComputation(self): # If we use LinBox, keep everything as-is. def sample(zeros): # Not currently sure how to handle this... maybe we'll just "reject" # for now, come back, and sub something else in later. We shouldn't # be halting computation. For now, we should raise an exception that # the Chain catches, and warns the user by exiting with exit code # 1 or 2. try: return np.array(ReducedKernelSample( self.matrices.coboundary, zeros, 2*self.dimension, self.faces, self.field, self._DEBUG ), dtype=self.dtypes[0]) except Exception as e: raise NumericalInstabilityWarning(e) self.sample = sample def _initial(self): """ Computes an initial state for the model's Complex. Returns: A numpy `np.array` representing a vector of spin assignments. """ return self.RNG.integers( 0, high=self.field, size=self.faces, dtype=self.dtypes[0] ) def _proposal(self, time): """ Proposal scheme for generalized Swendsen-Wang evolution on the Potts model. Args: time (int): Step in the chain. Returns: A numpy array representing a vector of spin assignments. """ # Compute the probability of choosing any individual cube in the complex. T = self.temperature(time) p = 1-np.exp(T) assert 0 <= p <= 1 # Choose cubes to include. uniform = self.RNG.uniform(size=self.cells) include = np.nonzero(uniform < p)[0] boundary = self.complex.Boundary[self.dimension] q = self.field # Evaluate the current spin assignment (cochain). coefficients = self.spins[boundary] coefficients[:,1::2] = -coefficients[:,1::2]%q sums = coefficients.sum(axis=1)%q zeros = np.nonzero(sums == 0)[0] includedZeros = np.intersect1d(zeros, include) # Sample from the kernel of the coboundary matrix, and evaluate again spins = self.sample(includedZeros) # If we're debugging, check whether the sample's correct. if self._DEBUG: coefficients = spins[boundary[includedZeros]] coefficients[:,1::2] = -coefficients[:,1::2]%q sums = coefficients.sum(axis=1)%q assert sums.sum() == 0 # Create output vectors. occupied = np.zeros(self.cells, dtype=self.dtypes[1]) occupied[includedZeros] = 1 coefficients = self.spins[boundary] coefficients[:,1::2] = -coefficients[:,1::2]%q sums = coefficients.sum(axis=1)%q zeros = np.nonzero(sums == 0)[0] satisfied = np.zeros(self.cells, dtype=self.dtypes[2]) satisfied[zeros] = 1 return spins, occupied, satisfied def _assign(self, cocycle): """ Updates mappings from faces to spins and cubes to occupations. Args: cocycle (np.array): Cocycle on the subcomplex. Returns: None. """ self.spins = cocycle