
# coding: utf-8

# In[1]:


get_ipython().magic('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()

