Percolation threshold simulation using C++

Question

I am studying Algorithms by Princeton University. The first week assignment is to simulate percolation by using Java.
https://coursera.cs.princeton.edu/algs4/assignments/percolation/specification.php
Since I am learning C++, I'm using C++ to do the assignment. Unfortunately, there is no grader for C++. Therefore, I came to here to ask for code review.

Core algorithm (from course website)

The core algorithm is weighted quick union with path compression.
The WeightedQuickUnionUF class represents a union–find data type (also known as the disjoint-sets data type). It supports the classic union and find operations, along with a count operation that returns the total number of sets.
The union–find data type models a collection of sets containing n elements, with each element in exactly one set. The elements are named 0 through n–1.Initially, there are n sets, with each element in its own set.The canonical element of a set (also known as the root) is one distinguished element in the set. Here is a summary of the operations:

• root(p) returns the canonical element(root) of the set containing p. It compress the path by updating the parent of p to its grandparent.
• The connected operation returns the same value for two elements if and only if they are in the same set.
• WeightedUnion(p, q) merges the set containing element p with the set containing element q. That is, if p and q are in different sets, replace these two sets with a new set that is the union of the two.
• count() returns the number of sets.

Percolation

• Model a percolation system using an n-by-n grid of sites.
• Each site is either open or blocked.
• A full site is an open site that can be connected to an open site in the top row via a chain of neighboring (left, right, up, down) open sites.
• The system percolates if there is a full site in the bottom row. In other words, a system percolates if we fill all open sites connected to the top row and that process fills some open site on the bottom row.

Goal

Find the threshold that the system will percolate. That is, the percentage of open site in the system for it to be percolated

WeightedQuickUnionUF.h

#pragma once

#include <vector>

class WeightedQuickUnionUF{
private:
    std::vector<int> parent; //parent link(site indexed)
    std::vector<int> sz; //size of component for roots(site indexed) / no of element in tree
    int count_; //number of components/tree
    
public:
    WeightedQuickUnionUF(int n);
    int count()const; //return no of component
    bool connected(int p, int q); //check if their components' are equal
    void WeightedUnion(int p, int q); //merge trees tgt, smaller tree becomes child of larger tree

private:
    int root(int p); //return root of p
};

WeightedQuickUnionUF.cpp

#include "WeightedQuickUnionUF.h"

WeightedQuickUnionUF::WeightedQuickUnionUF(int n){
    count_ = n;
    parent.reserve(n);
    for (int i = 0; i < n; i++){
        parent.push_back(i);
    }
    sz.reserve(n);
    for (int i = 0; i < n; i++){
        sz.push_back(i);
    }  
}

int WeightedQuickUnionUF::count()const{
    return count_;
}

int WeightedQuickUnionUF::root(int p){
    while (p != parent[p]){ //chase parent pointer until it reaches root => parent = parent's root => parent = root
        parent[p] = parent[parent[p]]; //change root of p to its grandparent: path-compression
        p = parent[p]; //go to grandparent
    }
    return p;
}

bool WeightedQuickUnionUF::connected(int p, int q){
    return root(p) == root(q);
}

void WeightedQuickUnionUF::WeightedUnion(int p, int q){
    int rootP = root(p);
    int rootQ = root(q);

    if (rootP == rootQ) return;
    if (sz[rootP] < sz[rootQ]){
        parent[rootP] = rootQ;
        sz[rootQ] += sz[rootP]; 
    }
    else {
        parent[rootQ] = rootP; 
        sz[rootP] += sz[rootQ]; 
    }
    --count_;
}

gen_uf.h

#pragma once
#include <random>

int rand_uf(int range_from, int range_to); //generate random number

gen_uf_im.cpp

#include "gen_uf.h"
#include <random>

int rand_uf(int range_from, int range_to){
    static std::random_device r1; //static: ensure only one seed in runtime
    static std::default_random_engine reng(r1()); //seed the random engine
    std::uniform_int_distribution<> distr(range_from, range_to); //initialize the range for the distribution
    return distr(reng); //transform the random unsigned int generated by reng into an int
}

Percolation.h

#pragma once
#include <vector>
#include <iostream>

#include "WeightedQuickUnionUF.h"
#include "gen_uf.h"

class Percolation: private WeightedQuickUnionUF{
private:
     //n^2 site
    std::vector<int> grid;

    //size of grid, 5-by-5 grid has size 5
    int sz_grid; 

    const int m_open = 1;

    const int closed = 0;

    //top virtual site
    int top; 

    //bottom virtual site
    int bottom;

    int numberOfOpenSites_;

public:
    //create n-by-n grid, with all sites initially closed
    Percolation(int n); 

    //open the site (row, col) if it is not open already
    void open(int row, int col); 
    
    // is the site (row, col) open? 
    bool isOpen(int row, int col)const; 

    // is the site (row, col) Full?
    bool isFull(int row, int col); 

    // return the number of open site
    int numberOfOpenSites()const; 

    // does the system percolates?
    bool percolates(); 

    //Monte Carlo Simulation
    double testPercolateThreshold();
};

Percolation.cpp

#include "Percolaton.h"
#include "WeightedQuickUnionUF.h"
#include "gen_uf.h"

#include <vector>
#include <iostream>

Percolation::Percolation(int n): 
    WeightedQuickUnionUF(n * n + 2)
{
    sz_grid = n;
    numberOfOpenSites_ = 0;
    grid.reserve(n * n);
    for (int i = 0; i < n * n; i++){ //n^2 complexity
        grid.push_back(closed); //initially all sites are closed
    }

    top = n * n; //top virtual site
    bottom = n * n + 1; //bottom virtual site

    int bot_tmp_ix = n * (n - 1) + 0; //bottom edge starts from [n - 1][0]
    
    for (int i = 0; i < n; i++){
        WeightedUnion(i, top); //connect top edge and top
        WeightedUnion(bot_tmp_ix, bottom); //connect bottom edge and bottom
        bot_tmp_ix++;
    }
} 

void Percolation::open(int row, int col){
    //row:1 col:0 for size 5 = grid[5]
    // = size(row) + col = 5(1) + 0
    int ix = sz_grid * row + col;
    if (!isOpen(row, col)) grid[ix] = m_open;
    ++numberOfOpenSites_;

    //connect to up, left, right, down sites if they are opened
    if ((row - 1) >= 0){ //up
        if(isOpen(row - 1, col)){
            int up_ix = sz_grid * (row - 1) + col;
            WeightedUnion(ix, up_ix);
        }
    }

    if ((col - 1) >= 0){ //left
        if(isOpen(row, col - 1)){
            int left_ix = sz_grid * row + (col - 1);
            WeightedUnion(ix, left_ix);
        }
    }

    if ((col + 1) <= sz_grid - 1){ //right
        if(isOpen(row, col + 1)){
            int right_ix = sz_grid * row + (col + 1);
            WeightedUnion(ix, right_ix);
        }
    }

    if ((row + 1) <= sz_grid - 1){ //down
        if(isOpen(row + 1, col)){
            int down_ix = sz_grid * (row + 1) + col;
            WeightedUnion(ix, down_ix);
        }
    }
}

bool Percolation::isOpen(int row, int col)const{
    int ix = sz_grid * row + col;
    if (grid[ix] == m_open) return true;
    else return false;
}

bool Percolation::isFull(int row, int col){
    int ix = sz_grid * row + col;
    return connected(ix, top);
}

int Percolation::numberOfOpenSites()const{
    return numberOfOpenSites_;
}

bool Percolation::percolates(){
    return connected(top, bottom);
}

double Percolation::testPercolateThreshold(){
    while (!percolates()){
        int row = rand_uf(0, sz_grid - 1);
        int col = rand_uf(0, sz_grid - 1); 
        if (!isOpen(row, col)){
            open(row, col);
        }
    }
    double threshold = numberOfOpenSites() / static_cast<double>(sz_grid * sz_grid);
    return threshold;
}

PercolationStat.h

#pragma once
#include "Percolaton.h"

class PercolationStat{
private:
    double thresholdSum;
    int size;
    int trials;

public:
    //perform independent trials on a n-by-n grid
    PercolationStat(int sz, int times);
    
    //sample mean of percolation threshold
    double PercolationMean();
    
};

PercolationStat.cpp

#include "PercolationStat.h"
#include "Percolaton.h"

PercolationStat::PercolationStat(int sz, int times): size(sz), trials(times)
{
    thresholdSum = 0;
    for (int i = 0; i < times; i++){
        Percolation p(sz);
        thresholdSum += p.testPercolateThreshold();
    }
}

double PercolationStat::PercolationMean(){
    return thresholdSum / trials;
}

ufmain.cpp

#include <iostream>
#include "WeightedQuickUnionUF.h"
#include "Percolaton.h"
#include "gen_uf.h"
#include "PercolationStat.h"

using namespace std;

int main(){
    int size;
    cout << "Enter the size(N) of percolation grid(N-by-N):";
    cin >> size;

    cout << "Enter number of time of simulation:";
    int time;
    cin >> time;

    PercolationStat pStat(size, time);
    cout << "Mean: " << pStat.PercolationMean() << endl;
}

Answer 1

Use std::size_t for sizes, counts and indices

You are using int to store sizes, counts and indices, however an int might not be large enough to store all possible sizes and indices of data structures that fit into memory. The proper type to use is std::size_t. Make a habit of using this instead of int. Note that this is also the type that is returned by many standard library container functions like std::vector::size(), and if you were to use int yourself, you would get compiler warnings when trying to compare those with the results of the STL functions.

Consider merging parent and sz into one std::vector

parent and sz should always have the same length, and the elements have to be in the same order. While it is working like you have written, you can enforce this by creating a small struct and make a single std::vector out of that:

class WeightedQuickUnionUF {
    struct Node {
        std::size_t parent;
        std::size_t sz;
    };

    std::vector<Node> nodes;
    ...
};

Then you can initialize it like so:

WeightedQuickUnionUF::WeightedQuickUnionUF(std::size_t n) {
    count_ = n;
    nodes.reserve(n);

    for (std::size_t i = 0; i < n; ++i) {
        nodes.push_back({i, i});
    }
}

And instead of writing root(p) you write nodes(p).root.

Naming things

While most names are chosen quite well, there are some inconsistencies. You have the issue that some member variable names clash with corresponding member function names. To avoid these clashes, in some cases you add an underscore at the end of the variable, in others you prefix the variable with m_. For non-clashing variable names you don't use any prefix or postfix. I recommend that instead you have a uniform way of naming all your member variables; either use the _ postfix or the m_ prefix.

Another issue is that some member functions start with a capital, others don't. It doesn't matter which style you use, but use it consistently. If you go for starting with a lower case letter, then of course the constructors are excluded from this rule as they have to match the name of the class exactly.

Composition versus inheritance

Inheritance is good when you have is-a relationships, for example a dog is an animal. However, can you say that a percolation is a weighted quick union algorithm? I recommend you don't let Percolation inherit from WeightedQuickUnionUF, but rather just make the latter a private member variable.

Storing grid more efficiently

For grid, you only store whether each grid site is open or closed. Thus, a bool would suffice. A bool normally only takes one byte instead of the typical 4 bytes for an int, but std::vector<bool> is often specialized such that it takes only one bit per element. (This has its drawbacks, but it's perfect for what you want to do.)

Don't use classes unnecessarily

Coming from Java you might think that everything needs to be in a class. However, in C++ you can have stand-alone functions, and I recommend you use them if you don't need a class. For example, PercolationStat does all its work in the constructor, and then you have a single member function to get the results out of it. This can be written much simpler as a single function instead:

double percolationStat(std::size_t sz, std::size_t times) {
    double thresholdSum = 0;

    for (std::size_t i = 0; i < times; ++i) {
        ...
    }

    return thresholdSum / times;
}

Related to this, I would move testPercolateThreshold() out of it and make it a stand-alone function. It can look like:

double testPercolateThreshold(std::size_t sz) {
    Percolation p(sz);

    while (!p.percolates()) {
        ...
    }

    return p.numerOfOpenSites() / static_cast<double>(sz * sz);
}

This reduces the responsibilties of class Percolation.

Use '\n' instead of std::endl

Prefer using '\n' instead of std::endl; the latter is equivalent to the former, but also forces the output to be flushed, which is usually not necessary and might hurt performance.

Unnecessary if-statements

Instead of if (condition) return true; else return false, just write return condition.

In open(), you check if isOpen() is false. However, in testPercolateThreshold() you do the same before calling isOpen(). So one of those tests is redundant. It might be safer to have open() return immediately if isOpen() is true, then it can be called unconditionally without corrupting the state.