Context:
--------

I had asked a [related question](https://stackoverflow.com/questions/50074975/what-would-be-an-efficient-way-to-find-the-averaged-size-distribution-without-ha) 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](https://www.ocf.berkeley.edu/~fricke/projects/hoshenkopelman/hk.c) 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.*)

[![enter image description here](https://i.stack.imgur.com/oeXJdm.jpg)](https://i.stack.imgur.com/oeXJdm.jpg)

* 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](https://www.ocf.berkeley.edu/~fricke/projects/hoshenkopelman/hk.c).