# A first Course in Random Matrix Theory
# by Marc Potters and Jean-Philippe Bouchaud
# Cambridge University Press 2020
# Python code used to generate graphs in Part II

%pylab inline
import scipy.linalg as npl
import scipy.optimize as opt
import scipy.special as special
import scipy.integrate as integrate
from cycler import cycler
from sklearn.isotonic import IsotonicRegression

Populating the interactive namespace from numpy and matplotlib

#Figures 14.1 14.2a and 14.a2b need the TracyWidom package written by Yao-Yuan Mao, Rutgers University, available at
#https://protect-eu.mimecast.com/s/-JMpC32RRFB65Rghgdl-C?domain=gist.github.com
#download the file TracyWidom.py and put it in the same directory as this notebook
try:
    import TracyWidom as TW
    TW_Available=True
except:
    TW_Available=False

ColorGraphs=True
if(ColorGraphs):
    HistoColor='powderblue'
else:
    plt.rcParams['axes.prop_cycle'] = cycler(color=['#000000', '#777777','#bbbbbb'])
    HistoColor='#777777'
plt.rc('legend', fontsize=8)
plt.rc('font', family='serif', size=10)
plt.rc('lines', linewidth=1, color='black')

#uncomment if latex is installed
plt.rc('text', usetex=True)
plt.rc('font', family='serif', serif='New Times Roman', size=10)
plt.rc('text.latex', preamble=r'\usepackage{mathptmx}\usepackage{amsfonts}')


#figure sizes in inches
figwidth=9./2.54
figwidth2=6/2.54
figheight=6/2.54
figdpi=200

#helper functions
myseed=180569
def wigner(myN,sigma2=1.):
    mat=np.random.randn(myN, myN) #non-symmetric Gaussian matrix
    return sqrt(sigma2/(2.*myN))*(mat+mat.transpose()) #symmetrize and normalize
def myg(values,z):#returns the sample Stieltjes transform of data values evaluated at points z (must be 1d array)
    return np.average(1./(z[:,newaxis]-values[newaxis,:]),axis=1)

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
random.seed(myseed)
dy=np.random.randn(1000)
y=cumsum(dy)*1./sqrt(1000)
x=arange(0,1,.001)
plt.plot(x,y)
plt.xlabel('$t$')
plt.ylabel('$B_t$')
plt.xlim(0,1)
#plt.ylim(0,0.345)
plt.tight_layout(pad=0.1)
#plt.savefig('fig8_1.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
random.seed(myseed)
nstep=50
Tmax=2000
delta=1./nstep
allx=np.zeros(nstep*Tmax)
xx=0
vol=1./sqrt(nstep)
for i in range(nstep*Tmax):
    vpo2dt=-xx/2.*delta
    xx=np.random.randn()*vol+vpo2dt+xx
    allx[i]=xx
Tmaxplot=40
x=arange(0,Tmaxplot,delta)
plt.plot(x,allx[:(Tmaxplot*nstep)])
plt.xlabel('$t$')
plt.ylabel('$X_t$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig8_2a.pdf',format='pdf')
plt.show()

random.seed(myseed)
plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
plt.hist(allx,100,density=True,color=HistoColor) 
x=arange(-4,4.01,.01)
y=exp(-x*x/2.)/sqrt(2*np.pi)
plt.plot(x,y)
plt.xlim(-4,4)
plt.xlabel('$X$')
plt.ylabel('$P(X)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig8_2b.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
nstep=2500
myN=25
mat=wigner(myN,.25)
allvall=np.zeros((nstep,myN))
for i in range(nstep):
    wval= npl.eigvalsh(mat)
    allvall[i]=wval
    mat=mat+wigner(myN,1./nstep)
allvall=allvall.transpose()
for j in range(myN):
    plt.plot(arange(0,1,1./nstep),allvall[j])
plt.xlabel('$t$')
plt.ylabel('$\lambda_t$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig9_1.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
nstep=100
Tmax=10
myN=25
mat=wigner(myN,.1)
allvall=np.zeros((nstep*Tmax,myN))
for i in range(nstep*Tmax):
    wval= npl.eigvalsh(mat)
    allvall[i]=wval
    nextval=wval-(0.5/nstep)*(wval+pow(wval,3)) #previous value plus restoring force for DBM in a potential
    mat=diag(nextval)+wigner(myN,1./nstep)
allvall=allvall.transpose()
for j in range(myN):
    plt.plot(arange(0,Tmax,1./nstep),allvall[j])
plt.xlabel('$t$')
plt.ylabel('$\lambda_t$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig9_2a.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
nstep=200
Tmax=10
myN=400
mat=wigner(myN,0.6) #initialize DBM with wigner matrix
allvall=np.zeros((nstep*(Tmax-3),myN))
for i in range(nstep*Tmax):
    wval= npl.eigvalsh(mat)
    if(i>=nstep*3): #only record eigenvalues after time t>3 to let the distribution attain steady-state
        allvall[i-nstep*3]=wval
    nextval=wval-(0.5/nstep)*(wval+pow(wval,3)) #previous value plus restoring force for DBM in a potential
    mat=diag(nextval)+wigner(myN,1./nstep) 
plt.hist(allvall.flatten(),100,density=True,color=HistoColor) 
def rho(x,gamma=0):
    a2=(sqrt(1+12*gamma)-1)/(6*gamma)
    a=sqrt(a2)
    def q(x):
        return gamma*x**2+1+2*gamma*a2
    return q(x)*sqrt(np.maximum(4*a2-x**2,0))/(2*np.pi)
x=arange(-1.6,1.6,.01)
y=rho(x,1.)
plt.plot(x,y)
plt.xlim(-1.6,1.6)
plt.xlabel('$\lambda$')
plt.ylabel('$\\rho(\lambda)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig9_2b.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=np.arange(0,6,.5)
y=x*0.
plt.plot(x,y,marker='x',linestyle='none',color='C0')
x=np.arange(7,8)
y=x*0.
plt.plot(x,y,marker='x',linestyle='none',color='C1')
x=np.arange(0,10)
y=x*0.
plt.plot(x,y,color='black')
myy=arange(-.5,.501,.01)
myx=6+sqrt(abs(1-4*myy**2))
plt.plot(myx,myy,color='C1',linestyle='--')
plt.annotate(text='', xy=(6,-1), xytext=(6,1), arrowprops=dict(arrowstyle='<-'))
plt.annotate(text='$\Lambda$', xy=(6,0), xytext=(5.5,.2), arrowprops=dict(arrowstyle='->'))
plt.annotate(text=r'$\mathfrak{z}(t)$', xy=(7,0), xytext=(7.4,.2), arrowprops=dict(arrowstyle='->'))
plt.annotate(text='$\lambda_{\max}$', xy=(5.5,0), xytext=(4.75,.2), arrowprops=dict(arrowstyle='->'))
plt.text(6.2,1,'$\Lambda+{\\rm i}\infty$',horizontalalignment='center')
plt.text(6.2,-1.1,'$\Lambda-{\\rm i}\infty$',horizontalalignment='center')
plt.xlabel('Re\,$z$')
plt.ylabel('Im\,$z$')
plt.xlim(0,8)
plt.ylim(-1.15,1.1)
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom='off',      # ticks along the bottom edge are off
    top='off',         # ticks along the top edge are off
    labelbottom='off',
    left='off',
    labelleft='off') # labels along the bottom edge are off

plt.tight_layout(pad=0.1)
#plt.savefig('fig10_1.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
eig=array([-1.34193837, -0.63666888, -0.23054108,  0.9994115 ,  1.33925136])
x=arange(-4.,4.,.01)
y=myg(eig,x)
plt.plot(x,y)
plt.ylim(-8,8)
plt.xlim(-4,4)
plt.xlabel('$z$')
plt.ylabel('$g_n(z)$')
plt.axvspan(eig[4], 10, alpha=0.25, color='gray')
plt.axvspan(-10, eig[0], alpha=0.25, color='gray')
plt.tight_layout(pad=0.1)
#plt.savefig('fig10_2a.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
delta=0.01
x1=arange(-12,eig[0],delta)
x2=arange(eig[4]+delta*.5,12,delta)
x=np.concatenate((x2,x1))
y=myg(eig,x)
plt.plot(y,x,label="$\mathfrak{z}(g)$")
x=eig[4]+0*y
plt.plot(y,x,linestyle="--",label="$\lambda_{\\rm max}$")
x=eig[0]+0*y
plt.plot(y,x,linestyle=":",label="$\lambda_{\\rm min}$")
plt.xlim(-8,8)
plt.ylim(-8,8)
plt.legend()
plt.xlabel('$g$')
plt.ylabel('$\mathfrak{z}(g)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig10_2b.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
def gwig(z):
    return (z-z*sqrt(1-4./z/z))/2
x=arange(-8,8,.01)-1e-7j
y=gwig(x).real
y2=gwig(x).imag
plt.plot(x.real,y,label="Re\,$\mathfrak{g}(z)$")
plt.plot(x.real,y2,label="Im\,$\mathfrak{g}(z)$")
plt.legend()
plt.xlabel('$z$')
plt.ylabel('$\mathfrak{g}(z)$')
plt.xlim(-8,8)
plt.axvspan(2, 10, alpha=0.25, color='gray')
plt.axvspan(-10, -2, alpha=0.25, color='gray')
plt.tight_layout(pad=0.1)
#plt.savefig('fig10_3a.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
xmax=4
x=arange(-1,1,.01)
y=1./x+x
plt.plot(x,y,label="$\mathfrak{z}(g)$")
x=arange(1,xmax+1,1)
y=2+0*x
plt.plot(x,y,linestyle="--",label="$\lambda_\pm=\pm 2$")
x=arange(-xmax,-.9,1)
y=-2+0*x
plt.plot(x,y,linestyle="--",color='C1')
x=arange(-xmax,-1,.01)
y=1./x+x
plt.plot(x,y,linestyle=":",label="wrong",color='C2')
x=arange(1,xmax+1,.01)
y=1./x+x
plt.plot(x,y,linestyle=":",color='C2')

plt.xlim(-xmax,xmax)
plt.ylim(-8,8)
plt.legend()
plt.xlabel('$g$')
plt.ylabel('$\mathfrak{z}(g)$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig10_3b.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
def rho3(x):
    return x*(x*x-3)/(2*np.pi*sqrt(4-x*x))
def rho4(x):
    x2=x*x
    return (x2*x2-4*x2+2)/(2*np.pi*sqrt(4-x2))
x=arange(-1.99,1.995,.01)
y=rho3(x)
plt.plot(x,y,label='$\delta\\rho_3(\lambda)$')
y=rho4(x)
plt.plot(x,y,label='$\delta\\rho_4(\lambda)$')

plt.legend()
plt.xlabel('$\lambda$')
plt.ylabel('$\delta\\rho(\lambda)$')
plt.ylim(-.5,.5)
plt.tight_layout(pad=0.1)
#plt.savefig('fig12_3.pdf',format='pdf')
plt.show()

def optbplus(a=4,n=10):
    nm2s=(n-2.)*(n-2.)
    return (nm2s*a+sqrt(nm2s*(n*n*a*a-4*n+4.)))/(2*(n-1.)*(n-2.))
def sad3(b,n,a):
    lam=a-b
    lam2=lam+n*b
    z=lam+1./lam
    return a*n*z-(n-1)*lam*lam/2.-lam2*lam2/2.-(n-1)*log(z-lam)-log(z-lam2)
def sad_dif(a,n):
    return sad3(optbplus(a,n),n,a)-sad3(0,n,a)
def t_star(n):
    n_min=2*sqrt(max(n-1.,0))/n
    return opt.brentq(sad_dif,n_min+1e-7,1.1,args=(n,))
t_starV=vectorize(t_star)

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
datan=arange(2.01,10,.01)
y=t_starV(datan)
y2=(2*sqrt(datan-1.)/datan)
if (plt.rcParams['text.usetex']):
    plt.plot(datan,y,label='$t_{\\rm\scriptsize c}$')
    plt.plot(datan,y2,label='$t_{\\rm\scriptsize s}$')
    plt.ylabel('$t_{\\rm\scriptsize c}$')
else:
    plt.plot(datan,y,label='$t_{\\rm c}$')
    plt.plot(datan,y2,label='$t_{\\rm s}$')
    plt.ylabel('$t_{\\rm c}$')
plt.legend()
plt.xlim(2,10)
plt.xlabel('$n$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig13_1.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
def myS(x):
    return np.heaviside(1-x,1)*x**2/2.+np.heaviside(x-1,0)*(2*x-log(x)-1.5)
myx=arange(1e-7,5,.05)
myy1=myS(myx)
myy2=0.5*myx*myx
myy3=2*myx
plt.plot(myx,myy1,label='quenched')
plt.plot(myx,myy2,label='annealed')
plt.plot(myx,myy3,linestyle='--',label='upper bound')
plt.xlabel('$t$')
plt.ylabel('$H_{\mathbf{X}}(t)$')
plt.legend()
plt.tight_layout(pad=0.1)
#plt.savefig('fig13_2.pdf',format='pdf')
plt.show()

if(TW_Available):
    plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
    x=arange(-5,3,.01)
    t1=TW.TracyWidom(1)
    t2=TW.TracyWidom(2)
    t3=TW.TracyWidom(4)
    y=t3.pdf(x)
    plt.plot(x,y,label="$f_4(u)$")
    y=t2.pdf(x)
    plt.plot(x,y,label="$f_2(u)$")
    y=t1.pdf(x)
    plt.plot(x,y,label='$f_1(u)$')
    plt.legend()
    plt.xlabel('$u$')
    plt.ylabel('$f_\\beta(u)$')
    plt.xlim(-5,3)
    plt.ylim(0,0.59)
    plt.tight_layout(pad=0.1)
    #plt.savefig('fig14_1.pdf',format='pdf')
    plt.show()
else:
    print("This plot requires TracyWidom package")

if(TW_Available):
    plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
    x=arange(-8,5,.01)
    ai=special.airy(x)
    y=ai[1]**2-x*ai[0]**2
    plt.plot(x,y,label='$\Phi_2(u)$')
    t2=TW.TracyWidom(2)
    y=t2.pdf(x)
    plt.plot(x,y,linestyle='--',label='$f_2(u)$')
    x=arange(-8,0.005,.01)
    y=sqrt(-x)/np.pi
    plt.plot(x,y,linestyle=':',label='$\sqrt{-u}/\pi$')
    plt.legend()
    plt.xlim(-8,2)
    plt.ylim(0,1)
    plt.xlabel('$u$')
    plt.tight_layout(pad=0.1)
    #plt.savefig('fig14_2a.pdf',format='pdf')
    plt.show()
else:
    print("This plot requires TracyWidom package")

if(TW_Available):
    plt.figure(figsize=(figwidth2,figheight),dpi=figdpi)
    x=arange(-3,8,.01)
    ai=special.airy(x)
    y=ai[1]**2-x*ai[0]**2
    plt.plot(x,y,label='$\Phi_2(u)$')
    t2=TW.TracyWidom(2)
    y=t2.pdf(x)
    plt.plot(x,y,linestyle='--',label='$f_2(u)$')
    x=arange(0.001,8,.01)
    y=exp(-(4./3)*x**(3./2))*13/96./x/np.pi
    plt.plot(x,y,linestyle=':',label='$\exp(-4u^{3/2}/3)/(8\pi u)$')
    plt.legend()
    plt.xlim(-3,7)
    plt.yscale('log')
    plt.ylim(1e-14,1)
    plt.xlabel('$u$')
    plt.tight_layout(pad=0.1)
    #plt.savefig('fig14_2b.pdf',format='pdf')
    plt.show()
else:
    print("This plot requires TracyWidom package")

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
def eigval_th(x):
    y=np.maximum(x,1.)
    return y+1./y 
random.seed(myseed)
myN=200
Ndata=200
xmax=4.
x=zeros(Ndata)
eigval=zeros(Ndata)
overlap=zeros(Ndata)
for k in range(Ndata):
    mat=wigner(myN)
    xx=(xmax*k)/Ndata
    mat[0][0]=mat[0][0]+xx
    val,vect=npl.eigh(mat,eigvals_only=False,eigvals=[myN-1,myN-1])#largest eigenvalue (with vector) of matrix mat
    x[k]=xx
    eigval[k]=val[0]
    overlap[k]=vect[0][0]**2
plt.plot(x,eigval,marker='.',linestyle='none',label='simulation',color='C1')
y2=eigval_th(x)
plt.plot(x,y2,label=r'$\lambda_1=\tilde{a}+1/\tilde{a}$; $\tilde{a}=\max(a,1)$',color='C0')
plt.legend()
plt.xlabel('$a$')
plt.ylabel('$\lambda_1$')
plt.xlim(0,4)
plt.tight_layout(pad=0.1)
#plt.savefig('fig14_3.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
plt.plot(x,overlap,marker='.',linestyle='none',label='simulation',color='C1')
def overlap_th(x):
    y=np.maximum(x,1.)
    return 1-y**(-2)
y=overlap_th(x)
if (plt.rcParams['text.usetex']):
    plt.plot(x,y,label='$|\mathbf{v}_1^{\scriptscriptstyle T}\mathbf{u}|^2=\max(1-a^{-2},0)$',color='C0')
    plt.ylabel('$|\mathbf{v}_1^{\scriptscriptstyle T}\mathbf{u}|^2$')
else:
    plt.plot(x,y,label='$|\mathbf{v}_1^{T}\mathbf{u}|^2=\max(1-a^{-2},0)$',color='C0')
    plt.ylabel('$|\mathbf{v}_1^{T}\mathbf{u}|^2$')
plt.legend()
plt.xlabel('$a$')
plt.xlim(0,4)
plt.ylim(0,1)
plt.tight_layout(pad=0.1)
#plt.savefig('fig14_5.pdf',format='pdf')
plt.show()

plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
x=arange(0.25,3,.05)
y=x+1./x
plt.plot(x,y,label='$g+1/g$')
x=[1]
y=[2]
plt.plot(x,y,marker='.',linestyle='none',label='$g_{+}$')
plt.legend()
plt.ylabel('$\mathfrak{z}(g)$')
plt.xlabel('$g$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig14_4.pdf',format='pdf')
plt.show()

import os
if(os.path.isfile("tracy_5.1000")):
    plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
    t1=TW.TracyWidom(1)
    x=arange(-15,15,.01)
    y=t1.pdf(x)*100
    x=2+x*.01
    data=np.loadtxt("tracy_5.1000")
    plt.hist(data,200,density=True,color=HistoColor,label='$N=1000$') 
    plt.plot(x,y,color='C0',label='Tracy-Widom')
    plt.plot(x,y)
    plt.legend()
    plt.yscale('log')
    plt.xlim(1.5,6)
    plt.ylim(1e-4,40)
    plt.xlabel('$\lambda_{\max}$')
    plt.ylabel('$P(\lambda_{\max})$')
    plt.tight_layout(pad=0.1)
    #plt.savefig('fig14_6.pdf',format='pdf')
    plt.show()
else:
    print("This figure needs a data file not provided with this notebook")

random.seed(myseed)
dataq=np.array([])
datavect=np.array([])
datavect2=np.array([])
for q in arange(.01,2,.025):
    T=int(floor(sqrt(1e6/q)))
    N=int(floor(q*T))
    aaa=np.random.randn(N, T)
    wei=np.heaviside(abs(aaa[0])-1,1)
    mat=wei*aaa
    bigm=aaa@mat.transpose()/T
    val,vect=npl.eigh(bigm,eigvals_only=False,eigvals=[N-1,N-1])#computes the largest eigenvalue (with vector) of bigm
    dataq=np.append(dataq,N/(float(T)))
    datavect=np.append(datavect,vect[0])
    wei=1-1./pow(aaa[0],2)
    mat=wei*aaa
    bigm=aaa@mat.transpose()/T
    val,vect=npl.eigh(bigm,eigvals_only=False,eigvals=[N-1,N-1])#computes the largest eigenvalue (with vector) of bigm
    datavect2=np.append(datavect2,vect[0])
plt.figure(figsize=(figwidth,figheight),dpi=figdpi)
p=special.erfc(1/sqrt(2))
c1=exp(-1/2)*sqrt(2/np.pi)+p
myx=arange(0,1,.01)
plt.plot(dataq,datavect2*datavect2,marker='s',markersize=2,linestyle='none',label='$f(y)=1-1/y$',color='C1')
plt.plot(dataq,datavect*datavect,marker='.',linestyle='none',label='$f(y)=\Theta(y-1)$',color='C0')
plt.plot(myx,1-c1/pow(c1-p,2)*myx,label='linear approx.',color='C0')
plt.ylim(0,1.1)
plt.xlim(0,2)
plt.legend()
plt.xlabel(r'$q$')
plt.ylabel(r'$\varrho$')
plt.tight_layout(pad=0.1)
#plt.savefig('fig14_7.pdf',format='pdf')
plt.show()