import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from matplotlib import rc

#graphical parameters
############################################################################### 
plt.style.use('ggplot')
rc('font',**{'family':'serif','serif':['Computer Modern Roman'],
     'size' : '12'})
rc('text', usetex=True)
rc('lines', linewidth=2)
plt.rcParams['axes.facecolor']='w'
###############################################################################

np.random.seed(42)

fig, axes = plt.subplots(nrows = 1, ncols = 3, figsize = (12,4))
H=[0.2,0.5,0.9]
hh=0
for ax in axes:
    for l in range(3):
        N= 1001 #number of samples
        H2= 2*H[hh] #Hurst parameter
        sigma=np.zeros(N)
        sigma[0]=1
        #constructs covariance generator 
        for k in range(1,N):
            sigma[k]=0.5* (  (k+1)**H2  -2.0*(k)**H2 + (k-1)**H2     )
            
        # we now construct the covariance matrix generator via circular embedding
        c=np.hstack((sigma,sigma[-2::-1][:-1])) 
        # notice that we don't assemble the covariance complete matrix
        labda=fft(c) # gets eigenvalues
        eta=np.sqrt(labda/(2*N))
        #samples complex normal vectors
        Z=np.random.standard_normal(len(c))+1.j*np.random.standard_normal(len(c))
        Zeta=Z*eta
        #computes fft
        X2n=fft(Zeta)
        A=X2n[:N+1]
        #gets real part
        X=np.real(A)
        c=1/N
        #rescales
        X1=c**(H2/2)*np.cumsum(X)
        #plots
        
        ax.plot(np.linspace(0,1,len(X1)),X1)
    ax.set_title('H = '+str(H2/2))
    plt.tight_layout()

    hh+=1
plt.show()