In [1]:
%matplotlib notebook
import numpy as np
from numpy import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#fun = lambda X,Y: 1*X**2 + 2*Y**2 + 3*X*Y + 4*X + 5*Y
fun = lambda X,Y: (np.sin(X)+np.cos(Y))**9
#fun = lambda X, Y: X**2 + Y**2 + 12
In [2]:
def gradient(a, b):
#val = [0, 0]
#val[0] = 2*a
#val[1] = 2*b
q = 9*(np.sin(a)+np.cos(b))**8
val = [q, q]
val[0] *= np.cos(a)
val[1] *= -np.sin(b)
return val
gr = lambda a, b: gradient(a, b)
In [3]:
def evol_alg(fun, eps, a, b, r_param):
point_a = [a]
point_b = [b]
point_f = [fun(a,b)]
old_p = fun(a,b)
eps = eps*r_param
change = 2*eps
max_iter = 200000
min_iter = 200
i = 1
while (change>eps and i<max_iter) or i<min_iter:
n_a = a + r_param*random.randn()
n_b = b + r_param*random.randn()
if fun(n_a, n_b)<old_p:
change = ((old_p-fun(n_a, n_b))**2)**(1/2) / (((n_a-a)**2+(n_b-b)**2) ** (1/2))
a = n_a
b = n_b
point_f.append(fun(a,b))
point_a.append(a)
point_b.append(b)
old_p = fun(a,b)
elif fun(n_a, b)<old_p:
change = ((old_p-fun(n_a, b))**2)**(1/2) / (((n_a-a)**2) ** (1/2))
a = n_a
point_f.append(fun(a,b))
point_a.append(a)
point_b.append(b)
old_p = fun(a,b)
elif fun(a, n_b)<old_p:
change = ((old_p-fun(a, n_b))**2)**(1/2) / (((n_b-b)**2) ** (1/2))
b = n_b
point_f.append(fun(a,b))
point_a.append(a)
point_b.append(b)
old_p = fun(a,b)
i += 1
if not(i<max_iter):
print("Zakończono przez przekroczenie maksymalnej liczby iteracji.")
elif not(change>eps):
print("Zakończono przez uzyskanie zmiany mniejszej od zadanego błędu.")
return a, b, old_p, point_a, point_b, point_f
In [4]:
def gradient_descent(fun, g_fun, eps, a, b, r_param):
point_a = [a]
point_b = [b]
point_f = [fun(a,b)]
old_p = fun(a,b)
eps *= r_param
change = 2*eps
max_iter = 20000
min_iter = 100
i = 0
while (change>eps and i<max_iter) or i<min_iter:
delta = g_fun(a, b)
n_a = a - r_param*delta[0]
n_b = b - r_param*delta[1]
change = ((old_p-fun(n_a, n_b))**2)**(1/2)
#change =((old_p-fun(n_a, n_b))**2)**(1/2) / (((n_a-a)**2+(n_b-b)**2) ** (1/2))
a = n_a
b = n_b
point_f.append(fun(a,b))
point_a.append(a)
point_b.append(b)
old_p = fun(a,b)
i += 1
if not(i<max_iter):
print("Zakończono przez przekroczenie maksymalnej liczby iteracji.")
elif not(change>eps):
print("Zakończono przez uzyskanie zmiany mniejszej od zadanego błędu.")
return a, b, old_p, point_a, point_b, point_f
In [5]:
def m_draw(ax, fun, point_a, point_b, point_f, m_color):
x_knots = np.linspace(min(point_a)-1,max(point_a)+1,120)
y_knots = np.linspace(min(point_b)-1,max(point_b)+1,120)
[X,Y] = np.meshgrid(x_knots, y_knots)
Z = fun(X,Y)
ax.plot_surface(X, Y, Z, cmap=plt.cm.rainbow)
ax.plot(point_a, point_b, point_f, '-*', color=m_color)
return ax
In [6]:
def m_draw_multi(ax, fun, point_a, point_b, point_f, m_color):
r1 = [min(sum(point_a, []))-1, max(sum(point_a, []))+1]
r2 = [min(sum(point_b, []))-1, max(sum(point_b, []))+1]
x_knots = np.linspace(r1[0],r1[1],1000)
y_knots = np.linspace(r2[0],r2[1],1000)
[X,Y] = np.meshgrid(x_knots, y_knots)
Z = fun(X,Y)
ax.plot_surface(X, Y, Z, cmap=plt.cm.rainbow)
for m in range(len(point_a)):
ax.plot(point_a[m], point_b[m], point_f[m], '-p', color=m_color[m])
return ax
In [7]:
x0 = 20*random.rand()-10
y0 = 20*random.rand()-10
ax = Axes3D(plt.figure())
a, b, old_p, point_a, point_b, point_f = evol_alg(fun, 0.0001, x0, y0, 0.1)
print([a, b])
ax = m_draw(ax, fun, point_a, point_b, point_f, 'black')
plt.show()
In [8]:
ax = Axes3D(plt.figure())
a, b, old_p, point_a, point_b, point_f = gradient_descent(fun, gr, 0.001, x0, y0, 0.0001)
print([a, b])
ax = m_draw(ax, fun, point_a, point_b, point_f, 'black')
plt.show()
In [9]:
ax = Axes3D(plt.figure())
foo = ['red', 'blue', 'black', 'yellow']
#foo = [(1.0, i, i) for i in np.linspace(0.8,0.1, 4)]
x = []
y = []
z = []
x0 = 20*random.rand()-10
y0 = 20*random.rand()-10
for m_color in foo:
a, b, old_p, point_a, point_b, point_f = evol_alg(fun, 0.0001, x0, y0, 0.1)
print([a, b])
x.append(point_a)
y.append(point_b)
z.append(point_f)
a, b, old_p, point_a, point_b, point_f = gradient_descent(fun, gr, 0.001, x0, y0, 0.0001)
print([a, b])
x.append(point_a)
y.append(point_b)
z.append(point_f)
foo.append('white')
ax = m_draw_multi(ax, fun, x, y, z, foo)
plt.show()
In [ ]: