Goal #5 of my dissertation proposal is to make ATEAMS++, my generalized model-sampling software, the industry standard. To me, that means:
- the algorithms are provably correct;
- it’s easy to use (including installation, troubleshooting, etc.);
- it has good documentation;
- it’s accompanied by a paper in a mathematical software/symbolic computation journal of resonable repute.
This paper takes care of #1, at least for the plaquette Swendsen–Wang (PSW) and plaquette invaded-cluster (PIC) algorithms. The remaining three items, however, have ensorcelled me for the past two weeks — so fascinated have I been that nearly all my working hours are spent on them. My therapist would tell me this is fine, but she’s on maternity leave and so is not in charge of me. I have decided that this is bad, and have come to believe that writing about it will exorcise whatever productivity demon hath possessed me.
Discrete probability is notoriously slow in its progress. It took nearly thirty years to prove that the critical probability of Bernoulli percolation on a square grid is \(\nicefrac 12 \). Sometimes, cutting through the murk — which more often than not means figuring out the right question to ask rather than answering an existing one — requires a computational touch.
The archetypal experiment in discrete probability goes like this:
- find out some qual/quantitative feature we want to measure or detect;
- come up with some code that simulates/samples the model on which we want to detect those features;
- run the sampler (or the Markov chain) for a long-ass time;
- compile the results.
It’s the scientific method. Pretty standard. There is not, however, standard software for doing these things — results in this vein are, most of the time, found using home-brewed programs that aren’t portable from system to system (or person to person). Since my first real research project as a PhD student was the above paper with Ben and Paul, I wanted to test my mettle and build a tool worth using. Unfortunately, I was (at the time) frightened of writing C++, and so tried to write everything in Python. This previous version of ATEAMS was extremely frustrating to install and to use. So, knowing that I could do far better with far less effort, I consigned the old ATEAMS to deprecation and embarked on writing ATEAMS++.
This time round, I want(ed, but I’m nowhere near done) to have a plan before monkeying around with code. The goal is to make ATEAMS++
- modular, so modules (e.g. arithmetic, different types of complexes, etc.) can be written and tested independently, and
- policy-based, so
- I have to write less code, and
- modules interface with each other infrequently and only in flexible, type-safe ways.
As a proponent of user stories, I started by sketching how a probabilist/statistical mechanist might structure an experiment on the plaquette random-cluster model at its self-dual point using the plaquette Swendsen–Wang algorithm (PSW). Mathematically, the \(d\)-dimensional plaquette random-cluster model on the CW complex \(X\) is a two-parameter family of probability measures \(\PRCM_{q,p}\) on the \(d\)-dimensional subcomplexes \(P \subset X\):
- \(q\) is the modulus of the ring \(\Z/q\Z\) in which our sampled (\(d-1\))-dimensional cochains take values, usually called the cluster weight;
- \(p\) is the Bernoulli trial density \(0 \leq p \leq 1\), i.e. the biased coin we flip to independently decide whether each compatible \(d\)-cell is included in \(P\); since we’re sampling at the self-dual point, we set \(p = \ln(\nicefrac{\sqrt q}{1 + \sqrt q})\).
(Note: as probability measures on subgraphs are called random graphs, probability measures on subcomplexes are called random complexes.)
Other models (like plaquette Bernoulli percolation) are parametrized in a similar way, suggesting a base programmatic workflow: users construct a complex, parametrize a model on that complex, then simulate the Markov chain underlying the model. The chain(model(complex)) nesting immediately tells us where the interface boundaries should lie — the complex knows nothing about the model, and the model knows nothing about the chain. It moreover says that both the complex and the model should be typed generically, so we get abstract Complex and Model classes from which actual complex/model definitions must inherit. The linear algebraic/topological subroutines can live in their own modules and be called on freely (even by the user), but primarily support the models.
A(n approximate) module diagram for ATEAMS++.
The snagged thread running through all the modules, however, was typing. Because each Complex represents a gigantic geometric object (e.g. a cubical complex) as a set of matrices and each Model calls on some support routine requiring those matrices (and because I didn’t want the code to be slow as shit) I wanted to use sparse linear algebra. Sadly most big-name sparse LA libraries don’t have support for finite rings/fields and the ones that do, in my uneducated view, are difficult to use. I thankfully found salvation in SparseRREF and PHAT, but this meant that I would have to do an untenable amount of method overloading, which wasn’t going to happen.
Enter the C++ template. Now, I could choose whether — and, importantly, how — the behavior of a class/method/function depends on the data it’s given. I also (perhaps mistakenly) took on a mathematician’s perspective and decided it would be better to abstractly define the ring over which the linear algebra is done. Each descendant of a generic Ring would be a parameter for each Model (and, regrettably, as a parameter for each Complex since the types need to be known at compile time). Here’s the templated definition:
typedef SparseRREF::field_t AbstractRing;
typedef unsigned long FINITE;
typedef SparseRREF::rat_t RATIONAL;
struct Ring {
public:
AbstractRing ring;
int characteristic;
};
(Note: before moving on, I want to stress that the code here is abridged from its published form.)
The first three lines just create more descriptive type aliases. AbstractRing is SparseRREF’s abstract field representation field_t that wraps custom data types (FINITE, for arithmetic over \(\Z/p\Z\)) with FLINT (RATIONAL, for arithmetic over \(\Q\)). Each Ring also has a characteristic property — for \(\Z/p\Z\) it’s \(p\), and for \(\Q\) it’s \(0\). I can then define some specific rings that I’ll use repeatedly, like \(\Z/p\Z\), \(\Q\), and \(\Z/2\Z\).
struct Q : public Ring {
public:
typedef RATIONAL dtype;
Q() {
this->ring = AbstractRing(SparseRREF::FIELD_QQ);
this->characteristic = 0;
};
Q(int characteristic) : Q() { };
};
struct Zp : public Ring {
public:
typedef FINITE dtype;
Zp(int characteristic) {
this->ring = AbstractRing(SparseRREF::FIELD_Fp, characteristic);
this->characteristic = characteristic;
};
};
struct Z2 : public Ring {
public:
typedef FINITE dtype;
Z2() {
this->ring = AbstractRing(SparseRREF::FIELD_Fp, 2);
this->characteristic = 2;
};
Z2(int characteristic) : Z2() { };
};
As I usually pass in the ring modulus/field characteristic as a command line argument, I overloaded the constructors for Q and Z2 to accept a superfluous parameter. This way, users (i.e. me) only have to change the initial declaration of the ring type, not the way the ring is instantiated. Now, as the support routines also need to know how each Complex wants its data stored, I added some SparseRREF-compatible matrix and vector template types.
typedef uint32_t INDEX;
template <typename RingLike>
using SparseMatrix = SparseRREF::sparse_mat<typename RingLike::dtype,INDEX>;
template <typename RingLike>
using SparseVector = SparseRREF::sparse_vec<typename RingLike::dtype,INDEX>;
template <typename RingLike>
using SparseMatrices = vector<SparseMatrix<RingLike>>;
template <typename T>
using DenseVector = std::vector<T>;
From here, we can define an abstract Complex class.
template <typename RingLike>
struct BoundaryData {
SparseMatrices<RingLike> Matrices;
SparseMatrix<RingLike> Full;
FlatBoundaryMatrix Flat;
};
template <typename RingLike>
class Complex {
public:
BoundaryData<RingLike> Boundary;
BoundaryData<RingLike> Coboundary;
bool periodic;
using RingType = RingLike;
virtual void constructBoundaryMatrices(Ring* R) = 0;
virtual void constructFlatBoundaryMatrix() = 0;
virtual void constructFullBoundaryMatrix(Ring* R) = 0;
virtual int size() = 0;
};
The boundary matrix-construction methods are left as virtual because they’re unique to each kind of complex. These methods are also made abstract (via the = 0 marker) to force any descendants to implement them. For example, the definition for a cubical complex looks like this:
template <typename RingLike>
class Cubical : public Complex<RingLike> {
public:
string name = "cubical";
vector<int> corners;
Cubical(vector<int> corners, bool periodic);
Cubical(vector<int> corners);
void constructBoundaryMatrices(Ring* R) override;
void constructFlatBoundaryMatrix() override;
void constructFullBoundaryMatrix(Ring* R) override;
int size() override;
};
Now, instantiating a complex representing the scale-\(10\) cubical \(4\)-torus \(\T^4_{10}\) with coefficients in the ring \(\Z/p\Z\) looks like
Cubical<Zp> plex({10,10,10,10});
That caps off the ATEAMS::complexes module. On to ATEAMS::models.
Because each Model is handed to some as-yet-undefined Markov chain simulator, we need ways to
- draw a sample from the model, and
- store the state of the model — i.e. the current time point and the sample just drawn.
Each Model thus needs a sample method that consumes, modifies, and returns a state-tracking object. And, because we don’t want the Model to interact up the hierarchy, the state-tracker should be owned by the chain, not the Model. The state-tracker should also be completely generic, in that it can be understood by every Model without specialization. That gives us something like
template <typename RingLike, template <typename> typename VectorLike>
struct ModelState {
VectorLike<RingLike> cochain;
vector<int> includes;
vector<int> essential;
int rank;
int energy;
int t;
};
Because different Model implementations store different kinds of state, the ModelState object acts as a catch-all; the only universally-modified member is t, the time point. That means the definition’s a bit clunky, but I got over it. Finally, as random complexes are (pretty much) ubiquitously parametrized by a (complex, coefficient ring, subcomplex dimension) triplet, the abstract Model class should look like
template <typename RingLike, template <typename> typename VectorLike>
class Model {
public:
Complex<RingLike>* complex;
Ring* coefficients;
int dimension;
bool DEBUG = false;
template <typename R>
using VectorType = VectorLike<R>;
using RingType = RingLike;
using State = ModelState<RingLike,VectorLike>;
Model(
Ring* R,
int dimension,
bool DEBUG=false
) : coefficients(R), dimension(dimension), DEBUG(DEBUG) { };
virtual State sample(int t, State& state) = 0;
virtual State initialize(State& state) = 0;
};
Abstracting the Ring and Complex classes was a good choice, since it lets the Model point to any objects — e.g. Zp and Cubical, respectively — that descend from them. I also added a DEBUG flag so I could figure out where I was messing up as I worked (and to help future ATEAMS++ developers do the same). Using a similar virtual ... = 0; pattern, every class that inherits from Model must implement a sample function that returns a ModelState of the appropriate type.
The plaquette Swendsen–Wang algorithm requires a bit more specificity than the base Model class provides, however. It needs to store both the \(p\) and \(q\) parameters. It also needs a source of (pseudo-)randomness that supports two distributions — one for the Bernoulli density \(p\) and another for sampling a new cochain. Since the \(q\) parameter is encoded by the ring R and I know the type of data the PSW algorithm stores, my SwendsenWang definition looks like this:
template <typename RingLike>
class SwendsenWang : public Model<RingLike,SparseVector> {
public:
double p;
string name = "Swendsen--Wang";
template <typename R>
using VectorType = SparseVector<R>;
using RingType = RingLike;
using State = ModelState<RingType,VectorType>;
SwendsenWang(
Complex<RingLike>* complex,
Ring* R,
int dimension,
double p;
bool DEBUG=false
);
State sample(int t, State& state) override;
State initialize(State& state) override;
private:
mt19937 RNG;
uniform_real_distribution<double> unituniform;
uniform_int_distribution<int> intuniform;
set<int> include;
};
Even though the plaquette Bernoulli percolation model is far simpler than PSW, its definition looks just about the same:
class Bernoulli : public Model<Z2,DenseVector> {
public:
double p;
string name = "Bernoulli";
template <typename R>
using VectorType = DenseVector<R>;
using RingType = Z2;
using State = ModelState<RingType,VectorType>;
Bernoulli(
Complex<Z2>* complex,
int dimension,
double p=0.5,
bool DEBUG=false
);
State sample(int t, State& state) override;
State initialize(State& state) override;
private:
std::mt19937 RNG;
std::uniform_real_distribution<double> unituniform;
std::vector<int> filtration;
};
Importantly, Model and its descendants expose scoped RingType, VectorType, and State type aliases that specify the return type of the sample method to the chain and the coefficient ring type to the support routines. Even better, the templates and aliases significantly clarify and shorten the code. All told, instantiating the SwendsenWang model with \(\Z/3\Z\) coefficients on the \(2\)-dimensional subcomplexes of \(\T^4_{10}\) looks like
using namespace ATEAMS;
using Complex = complexes::Cubical<Zp>;
using Model = models::SwendsenWang<Complex::RingType>;
Complex::RingType R(3);
Complex torus({10,10,10,10});
Model PSW(&torus, &R, 2, 0.63);
(Ignore the fact that it’s bad practice to use Model and Complex as type aliases.)
The last piece of the puzzle was the Markov chain simulator. Ripping off my collaborative work on GerryChain, I wanted to use a new C++23 feature: the generator. Much like a Python generator, a C++ generator spoofs a for loop over an object not necessarily in memory at runtime by repeatedly pausing and resuming a coroutine. Basically, instead of constructing an iterable container (like a vector) and then iterating over the container in sequence, a generator calls a special function to generate the next element of the sequence on the fly. This setup is perfect for simulating Markov chains since a Markov chain, by definition, only cares about what it sees right now, not what it’s seen before. It’s even better for simulating these Markov chains in particular — because the complexes are so big, trying to store all the samples would use up all the system memory.
Based on the Model and ModelState definitions, all that’s required for a Chain class is a way to store a pointer to the Model, the ModelState, and the total number of iterations. This is, in its entirety, the Chain class:
template <typename ModelType>
class Chain {
public:
using State = typename ModelType::State;
State state;
ModelType* model;
int steps;
Chain(ModelType* model, int steps) {
this->model = model;
this->steps = steps;
State state;
this->state = state;
};
generator<State> simulate() {
this->model->initialize(this->state);
for (int t=0; t < this->steps; t++) {
this->model->sample(t, this->state, this->options);
co_yield this->state;
}
}
};
The co_yield keyword tells the compiler that simulate() is a coroutine, which forms a simple wrapper function for managing state. That’s it. Running the PSW algorithm for 1,000 iterations on \(\T^4_{10}\) with \(\Z/3\Z\) coefficients looks like this:
using namespace ATEAMS;
using Complex = complexes::Cubical<Zp>;
using Model = models::SwendsenWang<Complex::RingType>;
using Chain = arithmetic::Chain<Model>;
Complex::RingType R(3);
Complex torus({10,10,10,10});
Model PSW(&torus, &R, 2, 0.63);
Chain M(&PSW, 1000);
for (Chain::State state : M.simulate()) {
<...>
}
Not too bad after all.
That’s basically everything in ATEAMS++. I tried to make the documentation as replete as possible, and hopefully it’s easy to install. Still working on the paper, though.