#!/usr/bin/env python
# coding: utf-8

# # Regressão Polinomial Com Descida de Gradiente
# 
# Nos notebooks [Regressão Linear Simples com Descida de Gradiente](https://nbviewer.jupyter.org/github/yurimalheiros/ai-notebooks/blob/master/ml/gdlinearregression.ipynb) e [Regressão Linear Múltipla com Descida de Gradiente](https://nbviewer.jupyter.org/github/yurimalheiros/ai-notebooks/blob/master/ml/gdmultiplelinearregression.ipynb) foi mostrado como ajustar funções lineares a um conjunto de dados para realizar previsões.
# 
# Em ambos os casos, uma função só se adequava bem aos dados de treinamento se eles seguissem uma tendência linear. 
# Entretanto, para dados que não seguiam essa tendência, a função não se adequava bem, o que acabaria gerando previsões ruins. Assim, para dados não lineares é possível ajustar outros tipos de funções que não sejam lineares através da Regressão Polinomial, gerando curvas ou superfícies mais complexas.

# Continuaremos utilizando o exemplo dos preços das casas de acordo com o seu tamanho. Abaixo temos o gráfico representando esses dados.

# In[1]:


from bokeh.io import show, output_notebook
from bokeh.plotting import figure
output_notebook()

sizes = [0.998, 1.500, 1.175, 1.232, 1.121, 0.988, 1.240, 1.501, 1.225, 1.552, 0.975, 1.121, 1.020,
        1.501, 1.664, 1.488, 1.376, 1.500, 1.256, 1.690, 1.820, 1.652, 1.777, 1.504, 1.831, 1.200 ]

prices = [25.9, 29.5, 27.9, 25.9, 29.9, 29.9, 30.9, 28.9, 35.9, 31.5, 31.0, 30.9, 30.0, 28.9, 36.9,
          41.9, 40.5, 43.9, 37.5, 37.9, 44.5, 37.9, 38.9, 36.9, 45.8, 41.0]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
show(p, notebook_handle=True)


# A Regressão Polinomial funciona de maneira muito parecida com a regressão linear, a principal diferença é o tipo da função que é gerada. A regressão linear gera retas ou superfícies representadas por funções do tipo: $f(x) = w_0x_0 + w_1x_1 + ... + x_2x_n$. Por outro lado, a Regressão Polinomial pode gerar curvas e superfícies dos mais variados tipos, por exemplo: $f(x) = w_0x_0 + w_1x_1^2$ ou $f(x) = w_0 + w_1x_1 + w_2x_2^3$.
# 
# Vamos fazer um exemplo tentando ajustar uma função do tipo: $f(x) = w_0 + w_1x^2$ aos dados de treinamento.
# 
# O código abaixo recebe os valores de x, $w_0$ e $w_1$ e calcula o resultado da função.

# In[2]:


def f(x, w0, w1):
    return w1*(x**2) + w0


# A função de erro médio quadrático continua a mesma. Ela é calculada através do código abaixo:

# In[3]:


def square_error(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)**2

        total += error

    return total/(2*len(X))


# Mais uma vez utilizaremos a descida de gradiente para encontrar a função que melhor se ajusta aos dados de treinamento. Com isso, durante o algoritmo, será necessário calcular as derivadas da função de erro médio quadrático.
# 
# O código abaixo calcula a derivada parcial em relação a $w_0$

# In[4]:


def d_error_w0(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)

        total += error

    return total/(len(X))


# E a função abaixo calcula a derivada parcial em relação a $w_1$

# In[5]:


def d_error_w1(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)

        total += error*x**2

    return total/(len(X))


# Com a função e as derivadas parciais em mãos, aplicaremos o algoritmo de descida de gradiente.
# No código abaixo, os pesos são atualizados 1000 vezes.

# In[6]:


w = [0, 0]
learning_rate = 0.1

for i in range(1000):        
    w0i = w[0] - learning_rate*d_error_w0(sizes, prices, w[0], w[1])
    w1i = w[1] - learning_rate*d_error_w1(sizes, prices, w[0], w[1])

    w[0] = w0i
    w[1] = w1i

print("Função gerada: {} + {}x^2".format(w[0], w[1]))


# No gráfico a seguir observarmos que a curva encontrada não é mais uma reta.

# In[11]:


import numpy as np

x_points = np.linspace(0.8, 2, 50)
y_points = [f(i, w[0], w[1]) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)


# Vamos fazer um segundo exemplo, agora tentando ajustar a função: $f(x) = w_3x^3 + w_2x^2 + w_1x + w_0$

# In[12]:


def f(x, w0, w1, w2, w3):
    return w3*(x**3) + w2*(x**2) + w1*x + w0


# Nesse exemplo, precisamos calcular quatro derivadas parciais, pois temos 4 coeficientes ($w_0, w_1, w_2, w_3$)

# In[13]:


def d_error_w0(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error

    return total/(len(X))

def d_error_w1(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x

    return total/(len(X))

def d_error_w2(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x**2

    return total/(len(X))

def d_error_w3(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x**3

    return total/(len(X))


# Por fim, realizando a descida de gradiente, encontramos os coeficientes.

# In[14]:


w = [0, 0, 0, 0]
learning_rate = 0.1

for i in range(1000):        
    w0i = w[0] - learning_rate*d_error_w0(sizes, prices, w[0], w[1], w[2], w[3])
    w1i = w[1] - learning_rate*d_error_w1(sizes, prices, w[0], w[1], w[2], w[3])
    w2i = w[2] - learning_rate*d_error_w2(sizes, prices, w[0], w[1], w[2], w[3])
    w3i = w[3] - learning_rate*d_error_w3(sizes, prices, w[0], w[1], w[2], w[3])

    w[0] = w0i
    w[1] = w1i
    w[2] = w2i
    w[3] = w3i
    
print("Função gerada: {} + {}x^2 + {}x^3 + {}x^4".format(w[0], w[1], w[2], w[3]))


# In[18]:


x_points = np.linspace(0, 3, 50)
y_points = [f(i, w[0], w[1], w[2], w[3]) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)