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 scipy import integrate
from math import pi,ceil,floor


def sqrt_g_star_function(k):
    #small useful function in GFF_Simulation
    New = (np.sin(k/2))**2
    return (np.sqrt((2/3)*np.sum(New,axis=1)))**(-1)



def GFF_Simulation(N):
    #We define the GFF using the method described in Sheffield "Gaussian free field for mathematicians", Section 4.4.
    vecs = np.mgrid[0:N, 0:N,0:N].reshape(3,-1).T    
    Lambda_Star = 2*(np.pi/N)*vecs
    Lambda_Star = Lambda_Star[1:] #(remove origin) 
    sqrt_g_star_k = sqrt_g_star_function(Lambda_Star)
    K = np.concatenate(([0], sqrt_g_star_k))

    Gaussians = np.random.normal(0,1, (N**3) )+ (1j)*np.random.normal(0,1, (N**3) )

    sdv_times_Gauss = K * Gaussians
    elements = sdv_times_Gauss.reshape((N,N,N))
    
    return np.real((N**(-1.5))* np.fft.fftn(elements))
        

def GFF_Image(N,h,percent,boundary,HTML):
    #This function generates an image of the clusters for the discrete GFF above level h, both as png and html
    #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='GFF'    
    pperco.Percolation_Image(N,GFF_Simulation(N)>h,h,percent,boundary,1,name,'h',HTML)



def GFF_Video(N,hrange,percent,boundary):
    #This function generates a video of the clusters for the discrete GFF, with frames corresponding to levels in hrange
    #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='GFF'
    pperco.Percolation_Video(N,GFF_Simulation(N),hrange,percent,boundary,1,name,'h')
       
    
def GFF_Slide(N,hrange,percent,boundary):
    #This function generates a slider of the clusters for the GFF above level h for all h in hrange
    #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='GFF'
    pperco.Percolation_Slide(N,GFF_Simulation(N),hrange,percent,boundary,1,name,'h')
    

def GFF_Density(h):
    #Computes the probability that the GFF is larger than h
    
    #This is a good approximation for the Green function at 0
    g=1.516386
    
    #First define the the density of the GFF at the origin
    f=lambda t : np.exp(-(t**2)/(2*g))/(np.sqrt(2*pi*g))
    return round(integrate.quad(f, h, np.inf)[0],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
hc=1.225 #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 h_c='+str(GFF_Density(hc)))
size_interval=0.525
hmin=hc-size_interval
hmax=hc+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(GFF_Density(hmax))+','+str(GFF_Density(hmin))+']')


#We first compute a few images as html and png. N needs to be small for the html to not be heavy
N=150
GFF_Image(N,hc,10,1,1)
hrange=np.arange(floor(hmin*10)/10,ceil(hmax*10)/10+0.1,0.1)
for h in hrange:
    GFF_Image(N,h,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
h_between_frames=round((hmax-hmin)/number_of_frames,5)
hrange=np.arange(hmin,hmax+h_between_frames/2,h_between_frames)
GFF_Video(N,hrange,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
hrange=np.arange(floor(hmin*10)/10,ceil(hmax*10)/10+0.1,0.1)
GFF_Slide(N,hrange,10,1)