# observed values of the inputs 𝐱 = (𝑥₁, …, 𝑥ᵣ) and outputs 
#                       outputs y = (y₁, …, yᵣ)
#
# Assumed that y values are "compliant" to the linear model (assumption):
#
#       y = 𝑓(𝐱) = 𝑏₀ + 𝑏₁𝑥₁ + ⋯ + 𝑏ᵣ𝑥ᵣ,
#
# Error function:
#
#      sum of squared residuals SSR = Σᵢ(𝑦ᵢ − 𝑓(𝐱ᵢ))² 
#
# Derived (normalized) cost function:
#
#      cost function 𝐶 = Σᵢ(𝑦ᵢ − 𝑏₀ − 𝑏₁𝑥ᵢ)² / (2𝑛).
#
# Gradient of C:
#
#   ∂𝐶/∂𝑏₀ = (1/𝑛) Σᵢ(𝑏₀ + 𝑏₁𝑥ᵢ − 𝑦ᵢ)    = mean(𝑏₀ + 𝑏₁𝑥ᵢ − 𝑦ᵢ)
#   ∂𝐶/∂𝑏₁ = (1/𝑛) Σᵢ(𝑏₀ + 𝑏₁𝑥ᵢ − 𝑦ᵢ) 𝑥ᵢ = mean((𝑏₀ + 𝑏₁𝑥ᵢ − 𝑦ᵢ) 𝑥ᵢ)


import numpy as np


# #############################################
# ssr_gradient(x, y, b):
#
#    x and y, which contain the observation inputs and outputs (np.array()
#    b that holds the current values of the decision variables 𝑏₀ and 𝑏₁. 
#
# Return: [best b[0] = mean(b[0] + b[1] * x[k] - y[k]), 
#          best b[1] = mean( (b[0] + b[1] * x[k] - y[k]) x[k] )]

def ssr_gradient(x, y, b):
    res = b[0] + b[1] * x - y		# Vector !!!
    return res.mean(), (res * x).mean() # .mean() is a method of np.ndarray

x = np.array([1, 2, 3])
y = np.array([1, 2, 3])
b = np.array([1,1])

print( ssr_gradient(x,y,b) )


def gradient_descent(gradient, x, y, start, learn_rate=0.1, n_iter=50, 
                     tolerance=1e-06):

    vector = start

    for _ in range(n_iter):
        diff = -learn_rate * np.array(gradient(x, y, vector))

        if np.all(np.abs(diff) <= tolerance):
            break
        vector += diff

    return vector


x = np.array([5, 15, 25, 35, 45, 55])
y = np.array([5, 20, 14, 32, 22, 38])

x = np.array([1, 2, 3])
y = np.array([1, 2, 3])

z = gradient_descent(ssr_gradient, x, y, start=[0.5, 0.5], learn_rate=0.01, 
		     n_iter=100_000)

print(z)