Percolation simulation: secretly stops at a certain recursion depth-CodePudding

I wrote a small code to do a phase transition simulation using percolation. The whole principle is explained in this video : https://www.youtube.com/watch?v=a-767WnbaCQ (btw one of the best math videos I ever saw). In order to understand the code, I'm afraid you will have to watch at least the first 5 minutes of the video first.

Here's a quick explanation of the code: The array main (which I guess is terribly named) is initialized as an empty 2d array. Every index of the array represents a node. The function group finds a "group" (a color, a phase) for every node (index in main) you enter. In order to do so, it checks if a link to another node is active, if the node is not part of the group, it adds it to the group and repeats the operation for this node. A group in the code is represented by a unique number in the main array. The links are represented by two arrays: connections_horizontal and connections_vertical.

The code works fine for p>0.5 but it starts doing something wonky for p<0.5. It's easy to see that the amount of recursions the code has to go through is enormous. If there's only one phase (one color), the recursion depth will be of size^2. Before I forced the code to change the max recursion depth I had a Stack Overflow error and the code crashed (at low p). Now, the code just stops for a low p (without any error message). My guess is that somehow python just stops at a certain recursion depth (at low p you have bigger phases so bigger recursion depth). Does anyone have an idea why? I apologize in advance if the code is confusing (I tried my best to comment it). Thank you in advance to anyone who has the courage to read this/watch the video.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import sys

np.random.seed(0)

size=200
x=size**2
sys.setrecursionlimit(x)
main=np.zeros((size,size)) #where everything happens, every index in main is a point on the plane
connections_horizontal=np.random.rand(size-1,size)  #horizontal connections between every points in main
connections_vertical=np.random.rand(size,size-1)  #vertical connections between every points in main
p=0.5 #parameter of percolation

def group(array,index,grp_nbr,p): #index of point for which we want to find the group, the group number (to distinguish the groups), parameter p (percolation)

    array[index]=grp_nbr
    # 1: find if the neighbouring connection is active
    # 2. if the connection is active, put the neighbouring point in the same group
    # 3. for every new point, repeat the procedure until either no connections, or every neighbour is in the same group

    if index[0]<size-1: #checking connections to the right
        if connections_horizontal[index]>p: #if connection is active
            if array[index[0] 1,index[1]]!=grp_nbr: #and neighbour is not in the same group
                array[index[0] 1,index[1]]=grp_nbr #adding the point to the group
                group(array,(index[0] 1,index[1]),grp_nbr,p) #repeat procedure for new point (recursion)
    if index[1]<size-1: #checking the connections up
        if connections_vertical[index]>p:
            if array[index[0],index[1] 1]!=grp_nbr:
                array[index[0],index[1] 1]=grp_nbr
                group(array,(index[0],index[1] 1),grp_nbr,p)
    if index[0]>0: #checking the connections to the left
        if connections_horizontal[index[0]-1,index[1]]>p:
            if array[index[0]-1,index[1]]!=grp_nbr:
                array[index[0]-1,index[1]]=grp_nbr
                group(array,(index[0]-1,index[1]),grp_nbr,p)
    if index[1]>0: #checking the connections down
        if connections_vertical[index[0],index[1]-1]>p:
            if array[index[0],index[1]-1]!=grp_nbr:
                array[index[0],index[1]-1]=grp_nbr
                group(array,(index[0],index[1]-1),grp_nbr,p)


for i in range(size): #finding a group for every point on the plane
    for j in range(size):
        if main[i,j]==0:
            group(main,(i,j),np.random.rand(1),p)

plt.imshow(main,cmap='rainbow',interpolation='none')
plt.savefig('test.png')
plt.show()

CodePudding user response：

The problem does not appear to be inherently with your algorithm, or the value of p specifically (although that allows for some improvement as well, but it is not the cause of the problem). Rather, it is with you playing with the recursion limit, which affects imshow.

Try running this:

import numpy as np
import matplotlib.pyplot as plt
import sys

sys.setrecursionlimit(7 ** 2)

main = np.array([[0.5666014542065752, 0.5666014542065752, 0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.5232480534666997, 0.09394051075844168], [0.5666014542065752, 0.5666014542065752, 0.2653894909394454, 0.2653894909394454, 0.5232480534666997, 0.5232480534666997, 0.5759464955561793], [0.5666014542065752, 0.5666014542065752, 0.2653894909394454, 0.2653894909394454, 0.5232480534666997, 0.5232480534666997, 0.5759464955561793], [0.2653894909394454, 0.9292961975762141, 0.2653894909394454, 0.2653894909394454, 0.5232480534666997, 0.5232480534666997, 0.5759464955561793], [0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.5232480534666997, 0.5232480534666997, 0.5759464955561793], [0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.31856895245132366, 0.5232480534666997, 0.6674103799636817], [0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.2653894909394454, 0.31856895245132366, 0.13179786240439217, 0.13179786240439217]])
plt.imshow(main,cmap='rainbow',interpolation='none')
plt.savefig('test.png')
plt.show()

That 7x7 array is generated if I run your script with:

size = 7
p = 0.5

And you'll note that it causes a recursion error even though none of your code is involved. The problem is caused within imshow. Even plt.imshow(main) will trigger the issue.

Try raising the recursion limit and the problem goes away:

sys.setrecursionlimit(52)

This is a value that happens to be sufficient for this example - but it's more than size ** 2, which works out to 49 for size = 7.

When imshow isn't causing trouble (which is likely the case for the values you shared in the example), the value of p isn't inherently a problem. For example, try this:

size=50
x=size**3
p=0.2

No problem at all - but increase the size to 60 and it fails (on my machine anyway). The problem is that your algorithm is very inefficient and for low p recursion depth can be exceeded or simply cause the system to run out of other resources. Recursion depth is limited for a reason: recursion is resource-intensive and raising the number doesn't magically get you more resources.

Anything you can write recursively can also be written iteratively - if you feel the algorithm you have is close to what you need, you should try to implement a non-recursive solution. But there's likely more efficient solutions that may work recursively.