Percolation in fortran

Introduction

I thought I should add another project¹ to my GitHub repositories and since I was mainly extending the Quantum Computing Simulator project I decided to implement something very simple, but still very rich, a model that despite being simple is still useful for many physical phenomena and also exhibits a phase transition… that is, something on Percolation theory. Because it is so simple, I spiced things a bit by changing the programming language, using Fortran instead of the usual C++. Why Fortran? Well, because Fortran is still a thing².

By the way, I added some things to the Quantum Computing project, including a Matrix Product State Simulator and a stabilizer formalism simulator able to execute Clifford gates… maybe I’ll add sometime some blog entries about those, too.

Theory

I don’t intend to describe it here, the Wikipedia page is a good start and also in the GitHub project’s README¹ I added a free book that contains much more than I could write here. Here is the book again: Percolation Theory Using Python³. It’s worth noting here that the code uses an algorithm called union-find for clusters identification.

Results

Here are some results for a 2×2 grid compared with the analytical results:

And here are the results for bigger and bigger grids:

Some code

The code on GitHub¹ is simple enough, one can look over it there, I don’t think it require much explanation.
I’ve decided that instead I should add here code in javascript that instead of computing statistics, it does something more eye-candy: it displays the clusters, colored.

So, here it is:

var perc = (function () {
    var canvas = document.getElementById("percolationCanvas");
    var ctx = canvas.getContext("2d");
    var percolationText = document.getElementById("percolationText");

    var percolation = {
        Size: 64,
        grid: [],
        probability: 0.5,
        clusters: [],
        displaySize: canvas.width / 64,

        Init: function () {
            var nr = this.Size * this.Size;
            this.probability = Math.random() / 8 + 0.5;
            this.grid = new Array(nr);
            for (var i = 0; i < nr; ++i) this.grid[i] = Math.random() < this.probability;
        },

        Find: function (parents, c) {
            while (parents[c] != c) {
                parents[c] = parents[parents[c]];
                c = parents[c];
            }
            return c;
        },

        UnionParents: function(parents, ranks, c1, c2) {
            c1 = this.Find(parents, c1);
            c2 = this.Find(parents, c2);
            if (c1 != c2) {
                if (ranks[c1] < ranks[c2])
                    [c1, c2] = [c2, c1];
                parents[c2] = c1;
                if (ranks[c1] == ranks[c2])
                    ++ranks[c1];
            }       
        },

        Union: function (parents, ranks, i1, j1, i2, j2, c) {
            var ind1 = i1 * this.Size + j1;
            var ind2 = i2 * this.Size + j2;
            if (i2 < 0 || i2 >= this.Size || j2 < 0 || j2 >= this.Size || !this.grid[ind1] || !this.grid[ind2])
                return c;

            var c1 = this.clusters[ind1];
            var c2 = this.clusters[ind2];
            if (c1 == 0 && c2 == 0) {
                this.clusters[ind1] = c;
                this.clusters[ind2] = c;
                parents[c] = c;
                ++c;
            }
            else if (c1 == 0 && c2 != 0)
                this.clusters[ind1] = c2;
            else if (c1 != 0 && c2 == 0)
                this.clusters[ind2] = c1;
            else
                this.UnionParents(parents, ranks, c1, c2);

            return c;
        },

        Compute: function () {
            var nr = this.Size * this.Size;

            this.clusters = new Array(nr).fill(0);
            var parents = new Array(nr).fill(0);
            var ranks = new Array(nr).fill(0);
            var c = 1;

            for (let i = 0; i < this.Size; ++i) {
                for (let j = 0; j < this.Size; ++j) {
                    if (this.grid[i * this.Size + j]) {
                        c = this.Union(parents, ranks, i, j, i, j + 1, c);
                        c = this.Union(parents, ranks, i, j, i + 1, j, c);
                    }
                }
            }

            for (let i = 0; i < this.Size; ++i) {
                for (let j = 0; j < this.Size; ++j) {
                    let ind = i * this.Size + j;
                    if (this.clusters[ind] != 0) this.clusters[ind] = this.Find(parents, this.clusters[ind]);
                    else if (this.grid[ind]) this.clusters[ind] = c++; // this is only for visual purposes, marks as clusters even clusters with a single site... does not exist in the fortran code
                }
            }
        },

        Percolates: function () {
            let offset = this.Size * (this.Size - 1);
            for (let i = 0; i < this.Size; ++i)
                for (let j = 0; j < this.Size; ++j)
                    if (this.clusters[i] != 0 && this.clusters[i] == this.clusters[offset + j])
                        return true;
            return false;
        },

        HSVtoRGB: function (h, s, v) {
            const C = v * s;
            const Hp = h / 60;
            const X = C * (1 - Math.abs(Hp % 2 - 1));
            let r = 0;
            let g = 0;
            let b = 0;

            if (Hp < 1) {
                r = C;
                g = X;
                b = 0;
            } else if (Hp < 2) {
                r = X;
                g = C;
                b = 0;
            } else if (Hp < 3) {
                r = 0;
                g = C;
                b = X;
            } else if (Hp < 4) {
                r = 0;
                g = X;
                b = C;
            } else if (Hp < 5) {
                r = X;
                g = 0;
                b = C;
            } else {
                r = C;
                g = 0;
                b = X;
            }

            const m = v - C;
            r += m;
            g += m;
            b += m;

            r = Math.floor(256 * r);
            g = Math.floor(256 * g);
            b = Math.floor(256 * b);

            return ((r << 16) | (g << 8) | b);
        },

        GetHue: function (r, g, b) {
            const V = Math.max(r, g, b);
            const C = V - Math.min(r, g, b);
            let h = 0;
            if (C != 0) {
                if (V == r) h = ((g - b) / C) % 6;
                else if (V == g) h = (b - r) / C + 2;
                else h = (r - g) / C + 4;
            }

            return 60 * h;
        },

        ClustersColors: function (clustersSet) {
            const colorsClusters = new Map();
            colorsClusters.set(0, '#FFFFFF');
            for (const cluster of clustersSet) {
                if (cluster == 0) continue;
                var unique = false;
                let cnt = 0;
                do {
                    unique = true;
                    const hue = 360 * Math.random();
                    const s = 0.5 + Math.random() / 2;
                    const v = 0.8 + Math.random() / 5;
                    let c = this.HSVtoRGB(hue, s, v);

                    colorsClusters.set(cluster, '#' + ('000000' + c.toString(16)).slice(-6));

                    for (const otherCluster of colorsClusters.keys()) {
                        if (otherCluster == cluster || otherCluster == 0) continue;
                        const ocolor = parseInt(colorsClusters.get(otherCluster).substr(1), 16);
                        const or = (ocolor >> 16);
                        const og = ((ocolor >> 8) & 0xFF);
                        const ob = (ocolor & 0xFF);
                        const ohue = this.GetHue(or, og, ob);
                        if (Math.abs(ohue - hue) < 1) {
                            unique = false;
                            break;
                        }
                    }

                    if (++cnt > 1000) break;
                } while (!unique);
            }
            return colorsClusters;
        },

        Display: function (ctx) {
            const clustersSet = new Set();
            for (let i = 0; i < this.Size; ++i)
                for (let j = 0; j < this.Size; ++j)
                    clustersSet.add(this.clusters[i * this.Size + j]);

            const colorsClusters = this.ClustersColors(clustersSet);

            let sz = this.displaySize * this.Size;
            ctx.clearRect(0, 0, sz, sz);
            for (let i = 0; i < this.Size; ++i)
                for (let j = 0; j < this.Size; ++j) {
                    let ind = j * this.Size + i;
                    ctx.fillStyle = colorsClusters.get(this.clusters[ind]);
                    ctx.fillRect(i * this.displaySize, j * this.displaySize, this.displaySize, this.displaySize);
                }

            // compute the value for 'probability' of what is actually in the grid
            let cnt = 0;
            for (let i = 0; i < this.Size; ++i)
                for (let j = 0; j < this.Size; ++j)
                    if (this.grid[i * this.Size + j]) ++cnt;

            let displayProb = cnt / (this.Size * this.Size);

            percolationText.innerHTML = 'Probability: ' + (displayProb * 100).toPrecision(3) + '%     ' + 'Percolates: ' + (this.Percolates() ? 'yes' : 'no');
        }
    };

    return function () {
        percolation.Init();
        percolation.Compute();
        percolation.Display(ctx);
    }
})();

perc();
setInterval(perc, 5000);

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

var perc = (function () {

var canvas = document.getElementById("percolationCanvas");

var ctx = canvas.getContext("2d");

var percolationText = document.getElementById("percolationText");

var percolation = {

Size: 64,

grid: [],

probability: 0.5,

clusters: [],

displaySize: canvas.width / 64,

Init: function () {

var nr = this.Size * this.Size;

this.probability = Math.random() / 8 + 0.5;

this.grid = new Array(nr);

for (var i = 0; i < nr; ++i) this.grid[i] = Math.random() < this.probability;

Find: function (parents, c) {

while (parents[c] != c) {

parents[c] = parents[parents[c]];

c = parents[c];

}

return c;

UnionParents: function(parents, ranks, c1, c2) {

c1 = this.Find(parents, c1);

c2 = this.Find(parents, c2);

if (c1 != c2) {

if (ranks[c1] < ranks[c2])

[c1, c2] = [c2, c1];

parents[c2] = c1;

if (ranks[c1] == ranks[c2])

++ranks[c1];

}

Union: function (parents, ranks, i1, j1, i2, j2, c) {

var ind1 = i1 * this.Size + j1;

var ind2 = i2 * this.Size + j2;

if (i2 < 0 || i2 >= this.Size || j2 < 0 || j2 >= this.Size || !this.grid[ind1] || !this.grid[ind2])

return c;

var c1 = this.clusters[ind1];

var c2 = this.clusters[ind2];

if (c1 == 0 && c2 == 0) {

this.clusters[ind1] = c;

this.clusters[ind2] = c;

parents[c] = c;

++c;

}

else if (c1 == 0 && c2 != 0)

this.clusters[ind1] = c2;

else if (c1 != 0 && c2 == 0)

this.clusters[ind2] = c1;

else

this.UnionParents(parents, ranks, c1, c2);

return c;

Compute: function () {

var nr = this.Size * this.Size;

this.clusters = new Array(nr).fill(0);

var parents = new Array(nr).fill(0);

var ranks = new Array(nr).fill(0);

var c = 1;

for (let i = 0; i < this.Size; ++i) {

for (let j = 0; j < this.Size; ++j) {

if (this.grid[i * this.Size + j]) {

c = this.Union(parents, ranks, i, j, i, j + 1, c);

c = this.Union(parents, ranks, i, j, i + 1, j, c);

}

for (let i = 0; i < this.Size; ++i) {

for (let j = 0; j < this.Size; ++j) {

let ind = i * this.Size + j;

if (this.clusters[ind] != 0) this.clusters[ind] = this.Find(parents, this.clusters[ind]);

else if (this.grid[ind]) this.clusters[ind] = c++; // this is only for visual purposes, marks as clusters even clusters with a single site... does not exist in the fortran code

}

Percolates: function () {

let offset = this.Size * (this.Size - 1);

for (let i = 0; i < this.Size; ++i)

for (let j = 0; j < this.Size; ++j)

if (this.clusters[i] != 0 && this.clusters[i] == this.clusters[offset + j])

return true;

return false;

HSVtoRGB: function (h, s, v) {

const C = v * s;

const Hp = h / 60;

const X = C * (1 - Math.abs(Hp % 2 - 1));

let r = 0;

let g = 0;

let b = 0;

if (Hp < 1) {

r = C;

g = X;

b = 0;

} else if (Hp < 2) {

r = X;

g = C;

b = 0;

} else if (Hp < 3) {

r = 0;

g = C;

b = X;

} else if (Hp < 4) {

r = 0;

g = X;

b = C;

} else if (Hp < 5) {

r = X;

g = 0;

b = C;

} else {

r = C;

g = 0;

b = X;

}

const m = v - C;

r += m;

g += m;

b += m;

r = Math.floor(256 * r);

g = Math.floor(256 * g);

b = Math.floor(256 * b);

return ((r << 16) | (g << 8) | b);

GetHue: function (r, g, b) {

const V = Math.max(r, g, b);

const C = V - Math.min(r, g, b);

let h = 0;

if (C != 0) {

if (V == r) h = ((g - b) / C) % 6;

else if (V == g) h = (b - r) / C + 2;

else h = (r - g) / C + 4;

}

return 60 * h;

ClustersColors: function (clustersSet) {

const colorsClusters = new Map();

colorsClusters.set(0, '#FFFFFF');

for (const cluster of clustersSet) {

if (cluster == 0) continue;

var unique = false;

let cnt = 0;

do {

unique = true;

const hue = 360 * Math.random();

const s = 0.5 + Math.random() / 2;

const v = 0.8 + Math.random() / 5;

let c = this.HSVtoRGB(hue, s, v);

colorsClusters.set(cluster, '#' + ('000000' + c.toString(16)).slice(-6));

for (const otherCluster of colorsClusters.keys()) {

if (otherCluster == cluster || otherCluster == 0) continue;

const ocolor = parseInt(colorsClusters.get(otherCluster).substr(1), 16);

const or = (ocolor >> 16);

const og = ((ocolor >> 8) & 0xFF);

const ob = (ocolor & 0xFF);

const ohue = this.GetHue(or, og, ob);

if (Math.abs(ohue - hue) < 1) {

unique = false;

break;

}

if (++cnt > 1000) break;

} while (!unique);

}

return colorsClusters;

Display: function (ctx) {

const clustersSet = new Set();

for (let i = 0; i < this.Size; ++i)

for (let j = 0; j < this.Size; ++j)

clustersSet.add(this.clusters[i * this.Size + j]);

const colorsClusters = this.ClustersColors(clustersSet);

let sz = this.displaySize * this.Size;

ctx.clearRect(0, 0, sz, sz);

for (let i = 0; i < this.Size; ++i)

for (let j = 0; j < this.Size; ++j) {

let ind = j * this.Size + i;

ctx.fillStyle = colorsClusters.get(this.clusters[ind]);

ctx.fillRect(i * this.displaySize, j * this.displaySize, this.displaySize, this.displaySize);

}

// compute the value for 'probability' of what is actually in the grid

let cnt = 0;

for (let i = 0; i < this.Size; ++i)

for (let j = 0; j < this.Size; ++j)

if (this.grid[i * this.Size + j]) ++cnt;

let displayProb = cnt / (this.Size * this.Size);

percolationText.innerHTML = 'Probability: ' + (displayProb * 100).toPrecision(3) + '% ' + 'Percolates: ' + (this.Percolates() ? 'yes' : 'no');

}

};

return function () {

percolation.Init();

percolation.Compute();

percolation.Display(ctx);

}

})();

perc();

setInterval(perc, 5000);

It’s not that complex as it might appear, a lot of this code just takes care of displaying the colors in such a way that the clusters are distinguishable for the human eye (a not so trivial thing, that’s why HSV is used).

What the code does

Here is the above code running:

0

Note that the ‘probability’ field displays the post hoc computed value out of the actual values, instead of the value for the (pseudo)random number generator, those two values can differ.

Conclusion

Another project is added to the collection… there are over 20 projects on GitHub now, this one being among the simplest ones. I don’t know yet what the next project will be about… I’m open to suggestions, so leave a message if you think of something interesting, computational physics related.

CodeProject

Percolation The project repository on GitHub ↩ ↩ ↩
Fortran is still a thing ↩
Percolation Theory Using Python An open access book from Springer ↩

Introduction

Theory

Results

Some code

What the code does

0

Conclusion

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply

Introduction

Theory

Results

Some code

What the code does

0

Conclusion

Share this:

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply