Finding black cluster size distribution for square grid at different probabilities of black elements (averaged over a certain number of iterations)

Question

Context:

I had asked a related question on Stack Overflow. On receiving several helpful hints from the commentators there, I could successfully implement averaging of the size distribution over a certain number of iterations. This question is about the modified version of the code I had posted on Stack Overflow.

Brief Description:

Consider a \$L \times L\$ size matrix \$M\$, whose entries can be either 0 or 1. Each element is 1 with probability \$p\$ and 0 with probability \$1 - p\$. I will call the elements labelled 1 as black elements and elements labelled 0 as white elements. I'm trying to write a code which:

• Generates a random matrix with entries 0's and 1's. I need to input size of matrix \$L\$ and the probability \$p\$.

• Labels all the black elements belonging to the same cluster with the same number. I'm using the basic framework of the Hoshen Kopelman algorithm for this purpose. I made some additions of my own so that even two black elements connected diagonally along a vertex are also considered to belong to the same cluster.

(Define a cluster of black elements as a maximal connected component in the graph of cells with the value of 1, where edges connect cells whose rows and columns both differ by at most 1 (so up to eight neighbours for each cell). In other words if two black elements of the matrix share an edge or a vertex consider them as belonging to the same black cluster. That is, think of the matrix as a large square and the elements as small squares in the large square.)

• Within a loop that runs from \$p = 0\%\$ to \$p = 100\%\$, the number of black clusters of each size is calculated and printed (averaged over a certain number of iterations as input by the user).

Code:

#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <time.h>
#include <string.h>
int *labels;
int  n_labels = 0;

int uf_find(int x) {
    int y = x;
    while (labels[y] != y)
        y = labels[y];

    while (labels[x] != x) {
        int z = labels[x];
        labels[x] = y;
        x = z;
    }
    return y;
}

int uf_union(int x, int y) {
    return labels[uf_find(x)] = uf_find(y);
}

int uf_make_set(void) {
    labels[0] ++;
    assert(labels[0] < n_labels);
    labels[labels[0]] = labels[0];
    return labels[0];
}

void uf_initialize(int max_labels) {
    n_labels = max_labels;
    labels = (int*)calloc(sizeof(int), n_labels);
    labels[0] = 0;
}

void uf_done(void) {
    free(labels);
    labels = 0;
}

#define max(a,b) (a>b?a:b)
#define min(a,b) (a>b?b:a)

int hoshen_kopelman(int **matrix, int m, int n) {

    uf_initialize(m * n / 2);

    for (int i = 0; i<m; i++)
        for (int j = 0; j<n; j++)
            if (matrix[i][j]) {

                int up = (i == 0 ? 0 : matrix[i - 1][j]);
                int left = (j == 0 ? 0 : matrix[i][j - 1]);

                switch (!!up + !!left) {

                case 0:
                    matrix[i][j] = uf_make_set();
                    break;

                case 1:
                    matrix[i][j] = max(up, left);
                    break;

                case 2:
                    matrix[i][j] = uf_union(up, left);
                    break;
                }
                int north_west;
                if (i == 0 || j == 0)
                    north_west = 0;
                else
                    north_west = matrix[i - 1][j - 1];

                int north_east;
                if (i == 0 || j == (n - 1))
                    north_east = 0;
                else
                    north_east = matrix[i - 1][j + 1];

                if (!!north_west == 1)
                {
                    if (i != 0 && j != 0)
                    {
                        if (matrix[i][j - 1] == 0 && matrix[i - 1][j] == 0)
                        {
                            if (!!matrix[i][j] == 0)
                                matrix[i][j] = north_west;
                            else
                                uf_union(north_west, matrix[i][j]);
                        }

                    }

                }
                if (!!north_east == 1)
                {
                    if (i != 0 && j != n - 1)
                    {
                        if (matrix[i - 1][j] == 0 && matrix[i][j + 1] == 0)
                        {
                            if (!!matrix[i][j] == 0)
                                matrix[i][j] = north_east;
                            else
                                uf_union(north_east, matrix[i][j]);
                        }
                    }
                }
            }
    int *new_labels = (int*)calloc(sizeof(int), n_labels);

    for (int i = 0; i<m; i++)
        for (int j = 0; j<n; j++)
            if (matrix[i][j]) {
                int x = uf_find(matrix[i][j]);
                if (new_labels[x] == 0) {
                    new_labels[0]++;
                    new_labels[x] = new_labels[0];
                }
                matrix[i][j] = new_labels[x];
            }

    int total_clusters = new_labels[0];

    free(new_labels);
    uf_done();

    return total_clusters;
}

void check_labelling(int **matrix, int m, int n) {
    int N, S, E, W;
    for (int i = 0; i<m; i++)
        for (int j = 0; j<n; j++)
            if (matrix[i][j]) {
                N = (i == 0 ? 0 : matrix[i - 1][j]);
                S = (i == m - 1 ? 0 : matrix[i + 1][j]);
                E = (j == n - 1 ? 0 : matrix[i][j + 1]);
                W = (j == 0 ? 0 : matrix[i][j - 1]);

                assert(N == 0 || matrix[i][j] == N);
                assert(S == 0 || matrix[i][j] == S);
                assert(E == 0 || matrix[i][j] == E);
                assert(W == 0 || matrix[i][j] == W);
            }
}

int cluster_count(int **matrix, int size, int N)
{
    int i;
    int j;
    int count = 0;
    for (i = 0; i < size; i++)
    {
        for (j = 0; j < size; j++)
        {
            if (matrix[i][j] == N)
                count++;
        }
    }
    return count;
}

int main()
{
    srand((unsigned int)time(0));
    int p = 0;

    printf("Enter number of rows/columns: ");
    int size = 0;
    scanf("%d", &size);

    printf("\n");
    FILE *fp;

    printf("Enter number of averaging iterations: ");
    int iterations = 0;
    scanf("%d", &iterations);

    int* M = (int*)calloc(size*size+1, sizeof(int));

    for (int p = 0; p <= 100; p++)
    {
        printf("\nAt p = %d:", p);
        for (int iterate = 1; iterate <= iterations; iterate++)
        {
            char str[100];
            sprintf(str, "BlackSizeDistribution%03i.txt", p);
            fp = fopen(str, "w");
            int **matrix;
            matrix = (int**)calloc(10, sizeof(int*));
            int** matrix_new = (int**)calloc(10, sizeof(int*));
            matrix_new = (int **)realloc(matrix, sizeof(int*) * size);
            matrix = matrix_new;
            for (int i = 0; i < size; i++)
            {
                matrix[i] = (int *)calloc(size, sizeof(int));
                for (int j = 0; j < size; j++)
                {
                    int z = rand() % 100;
                    z = z + 1;
                    if (p == 0)
                        matrix[i][j] = 0;
                    if (z < p)
                        matrix[i][j] = 1;
                    else
                        matrix[i][j] = 0;
                }
            }

            hoshen_kopelman(matrix, size, size);
            int highest = matrix[0][0];
            for (int i = 0; i < size; i++)
                for (int j = 0; j < size; j++)
                    if (highest < matrix[i][j])
                        highest = matrix[i][j];
            int* counter = (int*)calloc(sizeof(int*), highest + 1);
            int high = 0;
            for (int k = 1; k <= highest; k++)
            {
                counter[k] = cluster_count(matrix, size, k);
                if (counter[k] > high)
                    high = counter[k];
            }
            int* size_distribution = (int*)calloc(sizeof(int*), high + 1);

            for (int y = 1; y <= high; y++)
            {
                int count = 0;
                for (int z = 1; z <= highest; z++)
                    if (counter[z] == y)
                        count++;
                size_distribution[y] = count;
            }
            check_labelling(matrix, size, size);
            for (int i = 0; i < size; i++)
                free(matrix[i]);
            free(matrix);
            for (int k = 1; k <= high; k++)
            {
                M[k] = M[k] + size_distribution[k];
            }
        }
        for (int k = 1; k <= size*size; k++)
        {
            double temp1 = (double)M[k];
            double temp2 = (double)iterations;
            double temp = temp1 / temp2;
            if (M[k] > 0)
            {
                printf("\n%d\t%lf", k, temp);
            }
        }
        printf("\n");
        for (int i = 0; i < size * size; i++)
        {
            M[i] = 0;
        }
    }
}

Problem:

My code is working perfectly fine for matrices of size around 100. Takes around 10 seconds to generate complete output. However,it seems to take way too long for matrices of size around 1000 or for larger sizes. Any idea how to improve the execution time for this code? Any other suggestions to improve it in general, are also welcome.

P.S: I'm sorry there are no comment lines in the code. I will try to write them soon and edit those in. However, as a rough guideline I mostly reused the code implementation here.

Answer 1

I can certainly see ways to chip off a bit of time. But the nature of this algorithm is that it does not scale well with size. So,I doubt that you can get massive boost. It would be definitely beneficial to do the following:

• Allocate most of the structures only once instead of 100 times. None of them change in size with p. Moreover your code overwrites them.

• the hoshen_kopelman function already returns a value that you are ignoring and I think you are calculating it as 'highest' again just after the call.

• if you modify hoshen_kopelman function to output the counter[] and the size_distribution data, you will save going through that matrix several times.

• the check_label function looks like a debug function, you should remove that in the optimized build. Speaking of that, try to do O3 optimization on the compiler.

Getting counter[] and size_distribution in the hoshen_kopelman function (assuming you have sufficient memory allocated somewhere and initialized to zero:

u32 max_size = 0;
for (int i = 0; i<m; i++)
    for (int j = 0; j<n; j++)
        if (matrix[i][j]) {
            int x = uf_find(matrix[i][j]);
            if (new_labels[x] == 0) {
                new_labels[0]++;
                new_labels[x] = new_labels[0];
            }
            matrix[i][j] = new_labels[x];
            /* getting counter */
            counter[new_labels[x]]++;
           /* max size */
            max_size = (max_size > counter[new_labels[x]]) ? max_size : (counter[new_labels[x]]+1);
        }

        //calloc size_distribution
            for (int z = 1; z <= new_labels[0]; z++)
                //calloc already inits to zero
                size_distribution[counter[z]]++;

Answer 2

First, a very interesting question.

This review covers the C code, and may not answer the performance question, but with the suggestions it might be easier to find what is not scaling well.

When one copies existing code and modifies it, it might be better to improve the style of the code.

Style issues:

• Usage of Global Variables.
• Functions that are too long and too complex.

Global Variables
There are two global variables used in this program, n_labels and labels. In this program n_labels and labels can be defined in the function hoshen_kopelman() and passed into the uf_ functions that user them.

It is generally better to pass variables in where they are needed than to use global variables.

A major problem with global variables is writing and maintaining correct code. It is very difficult to see where global variables are modified and how they are used in different functions and files.

As programs get larger and additional source files are added global variables created in one file may conflict with global variables previously created in another file. This conflict will show up when trying to link multiple files. In the C programming language the way to avoid this is to declare variables that need to be global to the file using the static keyword so that the variables are limited in scope to the current file.

int uf_find(int x, int *labels) {
    int y = x;
    while (labels[y] != y)
        y = labels[y];

    while (labels[x] != x) {
        int z = labels[x];
        labels[x] = y;
        x = z;
    }
    return y;
}

int uf_union(int x, int y, int *labels) {
    return labels[uf_find(x, labels)] = uf_find(y, labels);
}

int uf_make_set(int *labels, int n_labels) {
    labels[0] ++;
    assert(labels[0] < n_labels);
    labels[labels[0]] = labels[0];
    return labels[0];
}

void uf_initialize(int max_labels, int *labels) {
    labels = (int*)calloc(sizeof(int), max_labels);
    labels[0] = 0;
}

void uf_done(int *labels) {
    free(labels);
    labels = NULL;
}

Functions that are Too Long and Too Complex
There are several functions that are hard to read because they are too complex and too long. The main() function is one example. This function can be broken up into at multiple functions.

By breaking the main() function up we expose some possible additional issues.

• Lack of error checking
• Unused resources

*NOTE the function averagingIteration(int p, int iterate, int M, int size) is still too complex.

static int getUserInput(int *size, int* iterations)
{
    int status = EXIT_SUCCESS;
    printf("Enter number of rows/columns: ");
    scanf("%d", size);

    if (*size < 0)
    {
        fprintf(stderr, "The number of rows and columns can not be less than 0.");
        status = EXIT_FAILURE;
    }

    printf("\n");

    printf("Enter number of averaging iterations: ");
    scanf("%d", &iterations);

    if (*iterations < 0)
    {
        fprintf(stderr, "The number averaging iterations can not be less than 0.");
        status = EXIT_FAILURE;
    }

    return status;
}

static int averagingIteration(int p, int iterate, int *M, int size)
{
    int status = EXIT_SUCCESS;
    FILE *fp = NULL;

    char str[100];
    sprintf(str, "BlackSizeDistribution%03i.txt", p);
    fp = fopen(str, "w");
    if (!fp)
    {
        status = EXIT_FAILURE;
        fprintf(stderr, "Can't open file %s for output\n", str);
        return status;
    }

    int **matrix = (int**)calloc(10, sizeof(*matrix));
    if (!matrix) {
        fprintf(stderr, "calloc of matrix failed in averagingIteration p = %d, iterate = %d\n", p, iterate);
        return EXIT_FAILURE;
    }

    int** matrix_new = (int**)calloc(10, sizeof(int*));
    if (!matrix_new) {
        fprintf(stderr, "calloc of matrix_new failed in averagingIteration p = %d, iterate = %d\n", p, iterate);
        return EXIT_FAILURE;
    }

    matrix_new = (int **)realloc(matrix, sizeof(int*) * size);
    if (!matix_new) {
        fprintf(stderr, "calloc of matix_new failed in averagingIteration p = %d, iterate = %d\n", p, iterate);
        return EXIT_FAILURE;
    }

    matrix = matrix_new;

    for (int i = 0; i < size; i++)
    {
        matrix[i] = (int *)calloc(size, sizeof(int));
        for (int j = 0; j < size; j++)
        {
            int z = rand() % 100;
            z = z + 1;
            if (p == 0)
                matrix[i][j] = 0;
            if (z < p)
                matrix[i][j] = 1;
            else
                matrix[i][j] = 0;
        }
    }

    hoshen_kopelman(matrix, size, size);

    int highest = matrix[0][0];
    for (int i = 0; i < size; i++)
        for (int j = 0; j < size; j++)
            if (highest < matrix[i][j])
                highest = matrix[i][j];

    int* counter = (int*)calloc(sizeof(int*), highest + 1);
    int high = 0;
    for (int k = 1; k <= highest; k++)
    {
        counter[k] = cluster_count(matrix, size, k);
        if (counter[k] > high)
            high = counter[k];
    }
    int* size_distribution = (int*)calloc(sizeof(int*), high + 1);

    for (int y = 1; y <= high; y++)
    {
        int count = 0;
        for (int z = 1; z <= highest; z++)
            if (counter[z] == y)
                count++;

        size_distribution[y] = count;
    }

    check_labelling(matrix, size, size);
    for (int i = 0; i < size; i++)
        free(matrix[i]);

    free(matrix);

    for (int k = 1; k <= high; k++)
    {
        M[k] = M[k] + size_distribution[k];
    }

    return status;
}

static int loopOverP(int p, int size, int iterations)
{
    int status = EXIT_SUCCESS;

    int* M = (int*)calloc(size*size + 1, sizeof(*M));
    if (!M) {
        fprintf(stderr, "calloc of M failed in loopOverP p = %d\n", p);
        return EXIT_FAILURE;
    }

    for (int p = 0; p <= 100; p++)
    {
        printf("\nAt p = %d:", p);
        for (int iterate = 1; iterate <= iterations; iterate++)
        {
            status = averagingIteration(p, iterate, M, size);
        }
        for (int k = 1; k <= size*size; k++)
        {
            double temp1 = (double)M[k];
            double temp2 = (double)iterations;
            double temp = temp1 / temp2;
            if (M[k] > 0)
            {
                printf("\n%d\t%lf", k, temp);
            }
        }
        printf("\n");
        for (int i = 0; i < size * size; i++)
        {
            M[i] = 0;
        }
    }

    return status;
}

int main()
{
    int status = EXIT_SUCCESS;

    srand((unsigned int)time(0));
    int size;
    int iterations = 0;

    status = getUserInput(&size, &iterations);
    if (status != EXIT_SUCCESS)
    {
        return status;
    }

    status = loopOverP(size, iterations);
    return status;
}

Error Checking

It is always a good idea to check user input to see that the value are valid, a negative size may cause errors in calloc() and a negative iterative count will prevent the iterate loop from executing.

Certain functions should always be checked to see if they return valid values the functions in this case are calloc() and fopen(). The function calloc() may return a NULL pointer if there is not enough memory to allocate. The function fopen() may return a NULL pointer for several reasons including insufficient privlege to create the file.

Resource Usage
In the function main() in this program multiple files are opened for output. These files are never written to in the program. This is a possible bug. None of these files are ever closed. In addition each call to fopen() is slowing down the program because fopen() is making system calls to open the files.

Answer 3

Any idea how to improve the execution time for this code?

Create a white border and make the matrix a bit larger.

Rather than have special "am I near the edge" code like cases i == 0 and j == (n - 1), just make the matrix (n+2*B)*(n+2*B) and fill the border with "white".

This will not reduce the O(), but will speed along code with the removed special cases.

I think the B can be 1, or maybe it needs to be 2.

   #define B 1

   for (int i = B; i<(m-2*B); i++)
    for (int j = B; j<(n-2*B); j++)
        if (matrix[i][j]) {

            // int up = (i == 0 ? 0 : matrix[i - 1][j]);
            int up = matrix[i - 1][j];
            // int left = (j == 0 ? 0 : matrix[i][j - 1]);
            int left = matrix[i][j - 1]); 
            ...

The marginal increase in memory size is linear with n and the extra initialization time trivial.

---

Any other suggestions to improve it in general, are also welcome.

Alternative allocation idiom.

In C, I find sizing to the de-referenced variable easier to code, review and maintain than sizing to the type. No need to code the matching type - and maybe get it wrong and no need to update.

// labels = (int*)calloc(sizeof(int), n_labels);
labels = calloc(sizeof *lables, n_labels);