In [4]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import collections, colors, transforms
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import sympy
from matplotlib import animation, rc
from IPython.display import HTML
import matplotlib.mlab as mlab
from scipy import stats
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
15.1 (a) Rosenbrock Function¶
In [2]:
n = 2
x = np.linspace(-n,n,1000)
X,Y = np.meshgrid(x,x)
Error = (1-X)**2 + 100*(Y-X**2)**2
fig =plt.figure(figsize=(6, 6), dpi= 180, facecolor='w', edgecolor='k')
plot = plt.imshow(np.log(Error),cmap = 'RdYlBu')
ax = fig.gca()
fig.colorbar(plot, shrink=0.8)
plt.title('log of the Rosenbrock function')
Out[2]:
15.1 (b) Downhill Simplex¶
In [201]:
##### PERFORM DH SIMPLEX #####
a = np.zeros((3,2))
#a[0,:]= np.random.rand(2)*1000-500
a[0,:] = -np.ones(2)
a[1,:] = a[0,:]
a[2,:] = a[0,:]
a[1,0] -= 1
a[2,1] -= 1
x = a[:,0]
y = a[:,1]
[xx,yy,err] = dhsimplex(x,y,1000)
n = 2
x = np.linspace(-n,n,100)
X,Y = np.meshgrid(x,x)
Error = (1-X)**2 + 100*(Y-X**2)**2
##### PLOT #####
fig = plt.figure()
fig=plt.figure(figsize=(8, 6), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, np.log(Error), rstride=1,
cstride=1, alpha=1, cmap = 'RdYlBu',
linewidth=0, antialiased=False, shade=True)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
for i in range(0,xx.shape[0]):
ax.plot(xx[i,0:2],yy[i,0:2],np.log(err[i,0:2]), color = 'k', linewidth = 0.5, marker = '.', markersize = 2)
ax.plot(xx[i,1:3],yy[i,1:3],np.log(err[i,1:3]), color = 'k', linewidth = 0.5, marker = '.', markersize = 2)
ax.plot(xx[i,[2,0]],yy[i,[2,0]],np.log(err[i,[2,0]]), color = 'k', linewidth = 0.5, marker = '.', markersize = 2)
plt.axis('off')
ax.view_init(elev=80., azim=40)
plt.title('Downhill Simplex Method')
plt.show()
fig=plt.figure(figsize=(3, 3), dpi= 180, facecolor='w', edgecolor='k')
plt.plot(np.mean(err,axis=1),'k',linewidth = 1)
plt.title('mean simplex error')
plt.show()
15.2 Spin Glass Simulated Annealing¶
In [243]:
ind = np.random.randint(0,n,1)
tmp = S
S[ind] *= -1
print(energy(J,S))
print(energy(J,tmp))
In [482]:
###### SIMULATED ANNEALING #####
n=100
tstep = 5000
S = np.sign(np.random.rand(n)-0.5)
J = np.random.randn(n)
clr = [plt.cm.RdYlBu(x) for x in np.linspace(0,0.4,3)]
fig=plt.figure(figsize=(8, 5), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
j = 0
for alpha in {0.1,0.01,0.001}:
#for i in range(0,10):
Shist,ehist,p,Emin = simAnneal(n,tstep,alpha,S,J)
ax.plot(ehist,linewidth = 3,color = clr[j],label = r'$\alpha = %s$' % alpha)
j+=1
ax.plot(np.ones(tstep)*Emin,'k--',linewidth = 1,label ='E min')
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Error')
plt.xlabel('Time')
plt.title('Simulated Annealing')
Out[482]:
In [ ]:
In [477]:
###### GENETIC ALGORITHM #####
n=100
tstep = 150
J = np.random.randn(n)
clr = [plt.cm.Spectral(x) for x in np.linspace(0,0.4,4)]
fig=plt.figure(figsize=(8, 5), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
j = 0
for beta in {10,1,0.1,0.01}:
#for i in range(0,5):
en,Emin = genetic(J,n,tstep,beta)
ax.plot(en,linewidth = 3,color = clr[j],label = r'$\beta = %s$' % beta)
j+=1
ax.plot(np.ones(tstep)*Emin,'k--',linewidth = 1,label ='E min')
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Error')
plt.xlabel('Time')
plt.title('Genetic Algorithm')
Out[477]:
FUNCTIONS¶
In [5]:
def energy(J,S):
return(-np.sum(J*S*np.roll(S,-1)))
In [391]:
def erf(X,Y):
return((1-X)**2 + 100*(Y-X**2)**2)
In [ ]:
def dhsimplex(x,y,stop):
xhist = np.zeros((stop,3))
yhist = np.zeros((stop,3))
err = np.zeros((stop,3))
j = 0
while j < stop:
tempx = []
tempy = []
xhist[j,:] = x
yhist[j,:] = y
for jj in range(0,3):
err[j,jj] = erf(x[jj],y[jj])
f = erf(x,y)
a = np.argsort(-f) #sort by error
x = x[a]
y = y[a]
xu = np.mean(x[1:3])
yu = np.mean(y[1:3])
if erf(2*xu-x[0],2*yu-y[0]) < erf(x[2],y[2]): #if the reflection is the new best point
if erf(3*xu-2*x[0],3*yu-2*y[0]) < erf(2*xu-x[0],2*yu-y[0]): # if growing gives a better point
x[0] = 3*xu-2*x[0]
y[0] = 3*yu-2*y[0]
#print('reflected and grew')
j+=1
continue
else:
#keep reflection
x[0] = 2*xu-x[0]
y[0] = 2*yu-y[0]
#print('reflected')
j+=1
continue
elif erf(x[1],y[1]) > erf(2*xu-x[0],2*yu-y[0]) > erf(x[2],y[2]): # if reflection is middle
#keep reflection
x[0] = 2*xu-x[0]
y[0] = 2*yu-y[0]
#print('reflected')
j+=1
continue
else: #reflecting gave the worst point
# try reflecting and shrinking
if erf((3/2)*xu-(1/2)*x[0],(3/2)*yu-(1/2)*y[0]) < erf(x[1],y[1]): # if reflecting and shrinking gives a better point
#keep it
x[0] = (3/2)*xu-(1/2)*x[0]
y[0] = (3/2)*yu-(1/2)*y[0]
#print('reflected and shrunk')
j+=1
continue
elif erf(0.5*(xu+x[0]),0.5*(yu+y[0])) > erf(x[1],y[1]): # perhaps just shrinking is better
x[0] = 0.5*(xu+x[0])
y[0] = 0.5*(yu+y[0])
#print('shrunk')
j+=1
continue
else: #shrink it all
tempx = x[2]
tempy = x[2]
for ii in range(0,x.shape[0]-1):
x[ii] = 0.5*(x[ii]+tempx)
y[ii] = 0.5*(y[ii]+tempy)
#print('shrunk everything')
j+=1
continue
return(xhist,yhist,err)
In [ ]:
def simAnneal(n,tstep,alpha,S,J):
ehist = np.zeros(tstep)
dE = np.zeros(tstep)
p = np.zeros(tstep)
bhist = np.zeros(tstep)
Shist = np.zeros((tstep,n))
Emin = -np.sum(np.abs(J))
for i in range(0,tstep):
# get current energy
ehist[i] = energy(J,S)
# get current spins
Shist[i,:] = S
# get tempurature parameter
beta = alpha*i
bhist[i] = beta
#flip a spin
ind = np.random.randint(0,n,1)
S[ind] *= -1
# get new energy
dE[i] = ehist[i] - energy(J,S)
p[i] = np.exp(-beta)
if dE[i] > 0:
continue
elif np.random.rand(1) < np.exp(-beta):
continue
else:
S[ind] *= -1 #adjust
return(Shist,ehist,p,Emin)
In [452]:
def genetic(J,n,tstep,beta):
Emin = -np.sum(np.abs(J))
tmp = np.zeros(n)
#p = np.zeros(n)
en = np.zeros(tstep)
S = np.sign(np.random.rand(n,n)-0.5)
for ii in range(0,tstep):
#get fitness etc
for i in range(n):
tmp[i] = energy(J,S[:,i])
#tmp = tmp/np.sum(tmp)
p = np.exp(-beta*(tmp-Emin))
#normalize p
p = p/np.sum(p)
t = np.cumsum(p)
en[ii] = np.min(tmp)
#create new strings
snew = np.zeros((n,n))
for j in range(0,n):
ind1 = np.where(np.random.multinomial(1,p))[0][0]
ind2 = np.where(np.random.multinomial(1,p))[0][0]
o = np.random.randint(0,100)
snew[0:o,j] = S[0:o,ind1]
snew[o:n,j] = S[o:n,ind2]
#mutate strings
r = [np.random.randint(0,100)]
snew[r,j] = -snew[r,j]
S = snew
#plt.plot(en,linewidth = 3)
#plt.plot(np.ones(tstep)*Emin,'k-')
#title('genetic algorithm')
return(en,Emin)
In [ ]:
fig, ax = plt.subplots()
patches = []
N = 99
for i in range(0,N+1):
p = np.vstack((xx[i,:],yy[i,:])).T
polygon = Polygon(p, True)
patches.append(polygon)
pat = PatchCollection(patches,cmap = 'Spectral', alpha=0.8)
colors = err
pat.set_array(np.array(colors))
ax.add_collection(pat)
ax.autoscale_view()
plt.show()
In [11]:
i = 1
# get current energy
ehist[i] = energy(J,S)
# get current spins
Shist[i,:] = S
# get tempurature parameter
beta = alpha*i
bhist[i] = beta
#flip a spin
ind = np.random.randint(0,n,1)
temp = S
temp[ind] = -temp[ind]
# get new energy
bhist[i] = np.exp(-beta*(e2-e1))
dE = energy(J,S) - energy(J,temp)
In [ ]:
for alpha in {0.1,0.01,0.001}:
#for i in range(0,5):
Shist,ehist,p,Emin = simAnneal(n,tstep,alpha,S,J)
ax.plot(ehist,linewidth = 2,color = clr[j],label = r'$\alpha = %s$' % alpha)
j+=1
ax.plot(np.ones(tstep)*Emin,'k--',linewidth = 1,label ='theoretical minimum')
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Error')
plt.xlabel('Time')
plt.title('Simulated Annealing')
#plt.suptitle(r"$\alpha = 0.01$")
#plt.subplot(2, 1, 1)
#plt.plot(ehist,linewidth = 0.7)
#plt.plot(np.ones(tstep)*Emin,'r--',linewidth = 0.8)
#plt.title('energy function')
#plt.subplot(2, 1, 2)
#plt.plot(p,color = 'orange')
#plt.title('temperature')
#fig.tight_layout()
#fig.subplots_adjust(top=0.8)