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.
Classes
class Bernoulli (C, dimension=1, PHAT=True)-
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.
PHAT:bool=True- Uses fast PHAT routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time.
Expand source code
class Bernoulli(): _name = "Bernoulli" def __init__(self, C, dimension=1, PHAT=True): """ 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. PHAT (bool=True): Uses fast PHAT routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time. """ # 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(PHAT) # Seed the random number generator. self.RNG = np.random.default_rng() def _delegateComputation(self, PHAT): if PHAT: 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) else: # If we can't/don't want to use LinBox, use inbuilt methods. Persistencer = Persistence(2, self.complex.flattened, self.dimension) def persist(filtration): return Persistencer.TwistComputePercolationEvents(filtration) 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)-
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.
Expand source code
class Glauber(): _name = "Glauber" def __init__( self, C, dimension=1, field=2, temperature=constant(-0.6), initial=None ): """ 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. """ 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,
stop=<function InvadedCluster.<lambda>>,
LinBox=True,
maxTries=16,
parallel=False,
minBlockSize=32,
maxBlockSize=64,
cores=4)-
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.
stop:function- A function that returns the number of essential cycles found before sampling the next configuration.
LinBox:bool=True- Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time.
maxTries:int=16- The number of attempts LinBox makes to sample a nonzero
vector in the kernel of the coboundary matrix A. See more
discussion in the
SwendsenWangmodel. parallel:boolean- Should matrix computations be done in parallel? (Uses C/C++).
minBlockSize:int=32- If
parallelis truthy, this is the smallest number of columns processed in parallel. maxBlockSize:int=64- If
parallelis truthy, this is the largest number of columns processed in parallel. cores:int=4- Number of available CPUs/cores/threads on the machine.
Included below are performance statistics for various configurations of
InvadedCluster. Each configuration completed 100 iterations on Pangolin, a Dell Precision 5280 workstation with an 18-core@1.305GhZ Intel Xeon W-2295 CPU. Warns indicates the number ofNumericalInstabilityWarnings caught during the run, and Zeros the number of all-zeros vectors returned.Samples in \mathbb T^2_N with PHAT s/iteration Field $N$ maxTries$\mathbb Z/2\mathbb Z$ 724 8 2.155 16 2.255 32 2.161 1024 8 4.921 16 4.663 32 4.915 1448 8 10.527 16 10.442 32 10.658 Expand source code
class InvadedCluster(): _name = "InvadedCluster" def __init__( self, C, dimension=1, field=2, initial=None, stop=lambda: 1, LinBox=True, maxTries=16, parallel=False, minBlockSize=32, maxBlockSize=64, cores=4 ): """ 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. stop (function): A function that returns the number of essential cycles found before sampling the next configuration. LinBox (bool=True): Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time. maxTries (int=16): The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix \(A\). See more discussion in the `SwendsenWang` model. parallel (boolean): Should matrix computations be done in parallel? (Uses C/C++). minBlockSize (int=32): If `parallel` is truthy, this is the smallest number of columns processed in parallel. maxBlockSize (int=64): If `parallel` is truthy, this is the largest number of columns processed in parallel. cores (int=4): Number of available CPUs/cores/threads on the machine. Included below are performance statistics for various configurations of `InvadedCluster`. Each configuration completed 100 iterations on Pangolin, a Dell Precision 5280 workstation with an 18-core@1.305GhZ Intel Xeon W-2295 CPU. **Warns** indicates the number of `NumericalInstabilityWarning`s caught during the run, and **Zeros** the number of all-zeros vectors returned. </br> </br> <center> <strong> Samples in \(\mathbb T^2_N\) with PHAT </strong> </center> ..include:: ./tables/InvadedCluster.PHATComputePersistencePairs.2.html """ # Object access. self.complex = C self.dimension = dimension self.stop = stop self._returns = 3 self.field = field # 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 < 10000 and LinBox: 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) self.nullity = len(self.complex.Boundary[self.dimension]) # Delegates computation for persistence and cocycle sampling. self._delegateComputation(LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries) # 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, LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries): 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. if LinBox: def sample(zeros): try: return np.array(LanczosKernelSample( self.matrices.coboundary, zeros, 2*self.dimension, self.faces, self.field, maxTries ), 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 LinBox and self.field > 2: def persist(filtration): essential = ComputePercolationEvents( self.matrices.full, filtration, self.dimension, self.field, self.complex.breaks ) essential = np.array(list(essential)) essential.sort() return essential[(essential >= low) & (essential < high)] # If we're using LinBox and the field we're computing over *is* two, # use PHAT. elif LinBox and 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) if not LinBox: # If we can't/don't want to use LinBox, use inbuilt methods. Reducer = MatrixReduction(self.field, parallel, minBlockSize, maxBlockSize, cores) Persistencer = Persistence(self.field, self.complex.flattened, self.dimension) coboundary = np.zeros((self.cells, self.faces), dtype=FINT) rows = self.matrices.coboundary[::3] cols = self.matrices.coboundary[1::3] entries = (self.matrices.coboundary[2::3]%self.field).astype(FINT) coboundary[rows,cols] = entries def persist(filtration): return Persistencer.TwistComputePercolationEvents(filtration) def sample(zeros): return KernelSample(Reducer, coboundary.take(zeros, axis=0)).astype(FINT) 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. """ # 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] stop = self.stop() occupied = np.zeros((self.rank, self.nullity)) satisfied = np.zeros(self.nullity) for t in essential: occupiedIndices = shuffledIndices[:t-low] occupied[j,occupiedIndices] = 1 if (j+1) == 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,
LinBox=True,
maxTries=16,
parallel=False,
minBlockSize=32,
maxBlockSize=64,
cores=4)-
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.
LinBox:bool=True- Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time.
maxTries:int=16- The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix.
parallel:boolean- Should matrix computations be done in parallel? (Uses C/C++).
minBlockSize:int=32- If
parallelis truthy, this is the smallest number of columns processed in parallel. maxBlockSize:int=64- If
parallelis truthy, this is the largest number of columns processed in parallel. cores:int=4- Number of available CPUs/cores/threads on the machine.
Expand source code
class Nienhuis(): _name = "Nienhuis" def __init__( self, C, q1, q2, dimension=2, field=2, initial=None, LinBox=True, maxTries=16, parallel=False, minBlockSize=32, maxBlockSize=64, cores=4 ): """ 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. LinBox (bool=True): Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time. maxTries (int=16): The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix. parallel (boolean): Should matrix computations be done in parallel? (Uses C/C++). minBlockSize (int=32): If `parallel` is truthy, this is the smallest number of columns processed in parallel. maxBlockSize (int=64): If `parallel` is truthy, this is the largest number of columns processed in parallel. cores (int=4): Number of available CPUs/cores/threads on the machine. """ # 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 < 10000 and LinBox: 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(LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries) def _delegateComputation(self, LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries): # If we use LinBox, keep everything as-is. if LinBox: def sample(faceZeros, cellZeros): try: return np.array(SubLanczosKernelSample( self.matrices.coboundary, cellZeros, faceZeros, self.field, maxTries=maxTries ), dtype=FINT) except Exception as e: raise NumericalInstabilityWarning(e) # If we use inbuilt routines, we need to actually construct the matrix # form of the coboundary matrix, which is extremely sparse (and really big). # Makes things slow. else: coboundary = np.zeros((self.cells, self.faces), dtype=FINT) rows = self.matrices.coboundary[::3] cols = self.matrices.coboundary[1::3] entries = (self.matrices.coboundary[2::3]%self.field).astype(FINT) coboundary[rows,cols] = entries Reducer = MatrixReduction(self.field, parallel, minBlockSize, maxBlockSize, cores) def sample(faceZeros, cellZeros): subbasis = coboundary.take(cellZeros, axis=0).take(faceZeros, axis=1) return KernelSample(Reducer, subbasis).astype(FINT) 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,
LinBox=True,
parallel=False,
minBlockSize=32,
maxBlockSize=64,
cores=4,
maxTries=16)-
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.
LinBox:bool=True- Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time.
parallel:bool=False- Should matrix computations be done in parallel? (Uses C/C++).
minBlockSize:int=32- If
parallelis truthy, this is the smallest number of columns processed in parallel. maxBlockSize:int=64- If
parallelis truthy, this is the largest number of columns processed in parallel. cores:int=4- Number of available CPUs/cores/threads on the machine.
maxTries:int=16- The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix.
If the first vector sampled by the Lanczos algorithm is all zeros, we perform the following steps until a nonzero vector is found, or we exhaust the number of attempts:
- sample once from A with no preconditioner;
- sample up to twice from D_0A, where D_0 is a random diagonal matrix;
- sample up to twice from A^\top D_1 A, where D_1 is a random diagonal matrix;
- sample for the remainder of the attempts from D_2 A^\top D_3 A D_2, where D_2, D_3 are random diagonal matrices.
If we spend the entire budget of attempts, it is likely that the result returned is either the all-zeros vector.
Included below are performance statistics for various configurations of
SwendsenWang. Each configuration completed 100 iterations on Pangolin, a Dell Precision 5280 workstation with an 18-core@1.305GhZ Intel Xeon W-2295 CPU. Warns indicates the number ofNumericalInstabilityWarnings caught during the run, and Zeros the number of all-zeros vectors returned.Samples in \mathbb T^2_N s/iteration Warns Zeros Field $N$ maxTries$\mathbb Z/2\mathbb Z$ 724 8 0.483 0 2 16 0.886 0 1 32 1.701 0 0 1024 8 1.190 0 3 16 1.969 0 1 32 3.351 0 0 1448 8 2.432 0 3 16 4.349 0 2 32 5.320 0 0 $\mathbb Z/3\mathbb Z$ 724 8 0.306 0 0 16 0.376 0 0 32 0.372 0 0 1024 8 0.751 0 0 16 0.915 0 0 32 0.919 0 0 1448 8 1.376 0 0 16 1.657 0 0 32 2.898 0 0 $\mathbb Z/5\mathbb Z$ 724 8 0.246 0 0 16 0.249 0 0 32 0.237 0 0 1024 8 0.551 0 0 16 0.567 0 0 32 0.506 0 0 1448 8 1.332 0 0 16 1.357 0 0 32 1.537 0 0 $\mathbb Z/7\mathbb Z$ 724 8 0.254 0 0 16 0.245 0 0 32 0.246 0 0 1024 8 0.529 0 0 16 0.537 0 0 32 0.526 0 0 1448 8 1.319 0 0 16 1.301 0 0 32 1.395 0 0 Samples in \mathbb T^4_N s/iteration Warns Zeros Field $N$ maxTries$\mathbb Z/2\mathbb Z$ 20 8 0.176 11 11 16 0.187 6 6 32 0.189 3 3 30 8 1.085 12 12 16 1.103 6 6 32 1.160 3 3 40 8 4.264 12 12 16 4.360 5 5 32 4.499 2 2 $\mathbb Z/3\mathbb Z$ 20 8 0.126 12 0 16 0.179 6 6 32 0.125 2 0 30 8 1.088 12 12 16 1.053 4 4 32 1.051 2 2 40 8 3.934 11 11 16 3.815 4 4 32 3.961 2 2 $\mathbb Z/5\mathbb Z$ 20 8 0.128 9 0 16 0.134 6 0 32 0.136 3 0 30 8 0.939 12 0 16 0.867 6 0 32 0.913 2 0 40 8 3.506 11 0 16 3.654 6 0 32 3.331 2 0 $\mathbb Z/7\mathbb Z$ 20 8 0.149 12 0 16 0.143 5 0 32 0.135 2 0 30 8 0.984 12 0 16 0.833 4 0 32 0.969 2 0 40 8 3.855 12 0 16 3.552 5 0 32 3.579 2 0 Expand source code
class SwendsenWang(): _name = "SwendsenWang" def __init__( self, C, dimension=1, field=2, temperature=constant(-0.6), initial=None, LinBox=True, parallel=False, minBlockSize=32, maxBlockSize=64, cores=4, maxTries=16 ): 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. LinBox (bool=True): Uses fast LinBox routines instead of slow inbuilt ones. WARNING: using inbuilt methods may dramatically increase computation time. parallel (bool=False): Should matrix computations be done in parallel? (Uses C/C++). minBlockSize (int=32): If `parallel` is truthy, this is the smallest number of columns processed in parallel. maxBlockSize (int=64): If `parallel` is truthy, this is the largest number of columns processed in parallel. cores (int=4): Number of available CPUs/cores/threads on the machine. maxTries (int=16): The number of attempts LinBox makes to sample a nonzero vector in the kernel of the coboundary matrix. If the first vector sampled by the Lanczos algorithm is all zeros, we perform the following steps until a nonzero vector is found, or we exhaust the number of attempts: 1. sample once from \(A\) with no preconditioner; 2. sample up to twice from \(D_0A\), where \(D_0\) is a random diagonal matrix; 3. sample up to twice from \(A^\top D_1 A\), where \(D_1\) is a random diagonal matrix; 4. sample for the remainder of the attempts from \(D_2 A^\top D_3 A D_2\), where \(D_2, D_3\) are random diagonal matrices. If we spend the entire budget of attempts, it is likely that the result returned is either the all-zeros vector. Included below are performance statistics for various configurations of `SwendsenWang`. Each configuration completed 100 iterations on Pangolin, a Dell Precision 5280 workstation with an 18-core@1.305GhZ Intel Xeon W-2295 CPU. **Warns** indicates the number of `NumericalInstabilityWarning`s caught during the run, and **Zeros** the number of all-zeros vectors returned. </br> </br> <center> <strong> Samples in \(\mathbb T^2_N\) </strong> </center> .. include:: ./tables/SwendsenWang.LanczosKernelSample.2.html </br> </br> <center> <strong> Samples in \(\mathbb T^4_N\) </strong> </center> .. include:: ./tables/SwendsenWang.LanczosKernelSample.4.html """ # Object access. self.complex = C self.temperature = temperature self.dimension = dimension self._returns = 2 self.field = field # 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 and LinBox: 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(LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries) def _delegateComputation(self, LinBox, parallel, minBlockSize, maxBlockSize, cores, maxTries): # If we use LinBox, keep everything as-is. if LinBox: 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(LanczosKernelSample( self.matrices.coboundary, zeros, 2*self.dimension, self.faces, self.field, maxTries=maxTries ), dtype=FINT) except Exception as e: raise NumericalInstabilityWarning(e) # If we use inbuilt routines, we need to actually construct the matrix # form of the coboundary matrix, which is extremely sparse (and really big). # Makes things slow. else: coboundary = np.zeros((self.cells, self.faces), dtype=FINT) rows = self.matrices.coboundary[::3] cols = self.matrices.coboundary[1::3] entries = (self.matrices.coboundary[2::3]%self.field).astype(FINT) coboundary[rows,cols] = entries Reducer = MatrixReduction(self.field, parallel, minBlockSize, maxBlockSize, cores) def sample(zeros): return KernelSample(Reducer, coboundary.take(zeros, axis=0)).astype(FINT) 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, dtype=FINT, size=self.faces ) 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[include]] coefficients[:,1::2] = -coefficients[:,1::2]%q sums = coefficients.sum(axis=1)%q zeros = np.nonzero(sums == 0)[0] # Create output vectors. satisfied = np.zeros(self.cells, dtype=FINT) satisfied[zeros] = 1 # Sample from the kernel of the coboundary matrix. spins = self.sample(zeros) return spins, 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