import numpy as np import plotly.io as pio pio.renderers.default='browser' import os os.chdir(os.path.dirname(os.path.abspath(__file__))) import Plot_Perco as pperco from math import ceil,floor #The random walk on the torus closes sites. This is marked by the number 0. Otherwise, the rest of the torus is open, i.e. 1. def RW(N,u, start): #Make sure the starting coordinates lie between 0 and N-1 t = int(np.floor(u * (N) **2 )) #The time the random walk runs for/N. We split the iteration in N parts for better performances. torus = np.ones((N,N,N),dtype=int) # The tuple current represents the position of the tuple during the walk current = start torus[start[0],start[1],start[2]] = 0 # We compute N times the successive steps made by the walk during time t. for i in range(N): DX = np.random.choice(6,t)-1 DX[DX>1]=0 # steps in x dir DY = np.random.choice(4,t)-1 DY[DY>1] = 0 # steps in y dir DY = DY*(1-abs(DX)) DZ = 2*np.random.choice(2,t)-1 # steps in z dir DZ = DZ*(1-abs(DX))*(1-abs(DY)) X = (current[0] + DX.cumsum()) % N # x positions Y = (current[1] + DY.cumsum()) % N # y positions Z = (current[2] + DZ.cumsum()) % N # z positions torus[X, Y, Z] = 0 current = [X[t-1],Y[t-1],Z[t-1]] end=current t = int(np.floor(u * (N) **3 ) - N * np.floor(u * (N) **2 )) #The time for the last few runs missing if t>0: DX = np.random.choice(6,t)-1 # steps in x dir DX[DX>1]=0 DY = np.random.choice(4,t)-1 # steps in y dir DY[DY>1] = 0 DY = DY*(1-abs(DX)) DZ = 2*np.random.choice(2,t)-1 # steps in z dir DZ = DZ*(1-abs(DX))*(1-abs(DY)) X = (current[0] + DX.cumsum()) % N # x positions Y = (current[1] + DY.cumsum()) % N # y positions Z = (current[2] + DZ.cumsum()) % N # z positions torus[X, Y, Z] = 0 #We also record the last point visited by the random walk end=[X[t-1],Y[t-1],Z[t-1]] return [torus,end] def RW_Output(N,urange): #This function generates a variable RW_All_Levels, which takes real values for each vertex of the torus, and such that RW_All_Levels>=u corresponds to all the points which have not been hit by the RW at time u(N^3), for any u in urange #We compute the first steps of the RW. [RW_Sample,end] = RW(N,urange[0],np.random.choice(N,3)) RW_All_Levels=RW_Sample*urange[-1]*3/2 #And now the rest of the RW. for k in np.arange(len(urange)-1): step=urange[k+1]-urange[k] [RW_Sample,end] = RW(N,step,end) #We remove the points which have been hit for the first time at this step, by assigning a value smaller than the current u, but bigger than the last step. RW_All_Levels[(RW_All_Levels>=urange[k]) & (RW_Sample==0)]=urange[k+1]-step/2 return RW_All_Levels def RW_Image(N,u,percent,boundary,HTML): #This function generates an image of the clusters for the RW on the torus at time uN^3, with frames corresponding to u in urange #Only the clusters with volume larger than (percent)/100 times the volume of the largest cluster are displayed #If percent is negative, we display only the (-percent) largest clusters instead. #If boundary=0 then the clusters have periodic boundary condition, and otherwise we use free boundary condition. #Put HTML=1 if you want to use the html file, otherwise put html=0 for better performance. name='RW' pperco.Percolation_Image(N,RW_Output(N,[u]),u,percent,boundary,1,name,'u',HTML) def RW_Video(N,urange,percent,boundary): #This function generates a video of the clusters for the RW on the torus at time uN^3, with frames corresponding to u in urange #Only the clusters with volume larger than (percent)/100 times the volume of the largest cluster are displayed #If percent is negative, we display only the (-percent) largest clusters instead. #If boundary=0 then the clusters have periodic boundary condition, and otherwise we use free boundary condition. name='RW' pperco.Percolation_Video(N,RW_Output(N,urange),urange,percent,boundary,1,name,'u') def RW_Slide(N,urange,percent,boundary): #This function generates a slider of the clusters for the RW on the torus at time uN^3, with frames corresponding to u in urange #Only the clusters with volume larger than (percent)/100 times the volume of the largest cluster are displayed #If percent is negative, we display only the (-percent) largest clusters instead. #If boundary=0 then the clusters have periodic boundary condition, and otherwise we use free boundary condition. name='RW' pperco.Percolation_Slide(N,RW_Output(N,urange),urange,percent,boundary,1,name,'u') def RW_Density(u): #Computes the probability that the RW at time uN^3 does not visit the origin #This is a good approximation for the Green function at 0 g=1.516386 return round(np.exp(-u/g),4) #%% os.chdir(os.path.dirname(os.path.abspath(__file__))) #Examples of parameters which produce nice outputs. N can be lowered to decrease computation time, and the time to display the html files uc=3.15 #This is a good approximation of the critical parameter #We display the density of the open vertices at the critical parameter print('Density of open edges at u_c='+str(RW_Density(uc))) size_interval=1.1 umin=uc-size_interval umax=uc+size_interval #Maximal and minimal value we want to display #We also display the size of the interval we chose in term of density, in order to be able to choose similar ones for different models print('Interval we plot in terms of density of open vertices=['+str(RW_Density(umax))+','+str(RW_Density(umin))+']') #We first compute a few images as html and png. N needs to be small for the html to not be heavy N=150 RW_Image(N,uc,10,1,1) urange=np.arange(floor(umin*10)/10,ceil(umax*10)/10+0.1,0.2) for u in urange: RW_Image(N,u,10,1,1) #We now make a video of the phase transition. As there is no html we can make N larger, but not too large due to plotly limitations N=200 number_of_frames=250 u_between_frames=round((umax-umin)/number_of_frames,5) urange=np.arange(umin,umax+u_between_frames/2,u_between_frames) RW_Video(N,urange,10,1) #Finally, we plot a few points as a 3D interactive slider. N needs to be quite small for the total html file to not be too heavy N=100 urange=np.arange(floor(umin*10)/10,ceil(umax*10)/10+0.1,0.2) RW_Slide(N,urange,10,1)