df = pd.read_csv("food_trucks.txt",header=None)
df.head()

plt.scatter(df[0], df[1])

y = np.array(df[1])
x1 = np.array(df[0])
x0 = np.ones(x1.shape)
print(len(y))
print(len(x1))
print(len(x0))

A = np.matrix([x0, x1]).transpose()
y = np.expand_dims(y, axis=1)
print(A.shape)
print(y.shape)

w = np.linalg.inv(A.transpose()*A)*(A.transpose()*y)

yhat = w[0]*x0 + w[1]*x1

plt.scatter(x1, np.array(yhat))
plt.scatter(x1, y, c='red')

w

"""
Explanation: МНК для линейной регрессии
$\mathbf{w}=(A^TA)^{-1}(A^T\mathbf{y})$ 
Загрузите файл food_trucks.txt. В нём два столбца значений — количество жителей в городе и доход грузовика с уличной едой в этом городе.
End of explanation
"""


from sklearn.metrics import r2_score

r2_score(np.squeeze(np.array(y)), np.squeeze(np.array(yhat)))

"""
Explanation: ${\displaystyle R^{2}}$ (Коэффициент детерминации) — это доля дисперсии объясняемая моделью
End of explanation
"""


f = lambda x: x ** 3 - 2*x ** 2 + 2
df = lambda x: 3 * x ** 2 - 4 * x # производная
x = np.linspace(-1, 2.5, 1000)
plt.plot(x, f(x))
plt.xlim([-1, 2.5])
plt.ylim([0, 3])
plt.show()

"""
Explanation: Предположения МНК

Истинная модель $y$ действительно линейна
Случайность выборки
Полнота ранга $A$
Случайность ошибок
Гомоскедастичность ошибок (однородность дисперсии)
Нормальность ошибок

МНК-оценки  являются несмещёнными, состоятельными и наиболее эффективными оценками в классе всех линейных несмещённых оценок
Методы оптимизации
Как было показано на лекции, большинство методов машинного обучения сводятся к поиску параметров, которые минимизируют ошибку на тренировочной выборке:
$$
\min_{w} L(w; D)
$$
Здесь:
* $D$ — размеченная обучающая выборка, ${x_i, y_i}_{i=1}^N$
* $L$ — функция потерь
* $w$ — настраиваемые веса алгоритма
В более общем виде задачу можно записать так:
$$
\min_{x} f(x)
$$
Здесь:
* $x$ — вектор значений
* $f$ — функция, принимающая вектор в качестве аргумента и выдающая числовое значение.
На семинаре рассмотрим подробнее методы минимизации функции, которые рассматривались на лекции.
Градиентный спуск
Для оптимизации возьмем простую функцию $f(x) = x^3 - 2x^2 + 2$
End of explanation
"""


def optimize_and_plot_steps(learning_rate, x_new=2, compute_learning_rate=None):
    x_old = 0
    # x_new — точка старта
    eps = 0.0001
    x_list, y_list = [x_new], [f(x_new)] # инициализируем список координат и значений функций при итерации
    
    # спускаемся, пока разница между координатами не достигла требуемой точности
    i = 0
    while abs(x_new - x_old) > eps: 
        x_old = x_new
        # считаем направление спуска
        direction = -df(x_old)
        # обновляем значение темпа обучения, если нам задана функция для этого
        if compute_learning_rate is not None:
            learning_rate = compute_learning_rate(i, learning_rate)
        # делаем шаг
        x_new = x_old + learning_rate * direction
        # запоминаем очередной шаг минимизации
        x_list.append(x_new)
        y_list.append(f(x_new))
        i += 1
        
    print("Found local min:", x_new)
    print("Steps number:", len(x_list))
    
    plt.figure(figsize=[10,3])
    
    plt.subplot(1,2,1)
    plt.scatter(x_list, y_list, c="r")
    plt.plot(x_list, y_list, c="r")
    plt.plot(x, f(x), c="b")
    plt.xlim([-1,2.5])
    plt.ylim([0,3])
    plt.title("Descent trajectory")

    plt.subplot(1,2,2)
    plt.scatter(x_list,y_list,c="r")
    plt.plot(x_list,y_list,c="r")
    plt.plot(x,f(x), c="b")
    plt.xlim([1.2,2.1])
    plt.ylim([0,3])
    plt.title("Descent trajectory (zoomed in)")
    plt.show()

"""
Explanation: И определим функцию, которая будет оптимизировать функцию $f(x)$ градиентным спуском с заданным постоянным шагом (он же learning rate, темп обучения).
End of explanation
"""


optimize_and_plot_steps(0.1, x_new=1)

"""
Explanation: Попробуем оптимизацию с шагом 0.1
End of explanation
"""


optimize_and_plot_steps(0.4)

"""
Explanation: Возьмем шаг побольше.
End of explanation
"""


optimize_and_plot_steps(0.5)

"""
Explanation: Что, если взять 0.5?
End of explanation
"""


optimize_and_plot_steps(0.49)

"""
Explanation: Застопорились в нуле, т.к. нашли точный локальный максимум. В нем производная равна нулю и мы никуда не можем сдвинуться. А если взять 0.49?
End of explanation
"""


optimize_and_plot_steps(0.51)

"""
Explanation: Что, если взять 0.51?
End of explanation
"""


optimize_and_plot_steps(0.05)

"""
Explanation: Мы улетели далеко влево. Это можно понять, распечатав значения x_new.
Теперь возьмём маленький шаг. Например, 0.05.
End of explanation
"""


optimize_and_plot_steps(0.01)

"""
Explanation: 0.01?
End of explanation
"""


def compute_learning_rate(i, prev_lr):
    return prev_lr * 0.9

optimize_and_plot_steps(0.45, compute_learning_rate=compute_learning_rate)

"""
Explanation: Чем меньше шаг, тем медленнее мы идём к минимум (и можем вдобавок застрять по пути). Чем больше темп обучения, тем большие расстояния мы перепрыгиваем (и имеем гипотетическую возможность найти минимум получше). Хорошая стратегия — начинать с достаточно большого шага (чтобы хорошо попутешествовать по функции), а потом постепенно его уменьшать, чтобы стабилизировать процесс обучения в каком-то локальном минимуме.
Теперь будем изменять шаг динамически:
$lr(i + 1) = lr(i) * 0.9$.
End of explanation
"""


optimize_and_plot_steps(0.45)

"""
Explanation: Если сравнивать с постоянным темпом обучения, то мы нашли минимум в 2 раза быстрее.
End of explanation
"""


df = pd.read_csv("food_trucks.txt", header=None)

df.head()

plt.scatter(df[0], df[1])

"""
Explanation: Это, конечно, искуственный пример, но такая же идея используются для обучения алгоритмов машинного обучения с миллионами параметров, функции потерь которых имеют очень сложную структуру и не поддаются визуализации.
Настройка линейной регрессии с помощью градиентного спуска
Рассмотрим теперь реальные данные и попробуем использовать градиентный спуск для решения задачи линейной регрессии.
End of explanation
"""


y = np.array(df[1])
x1 = np.array(df[0])
x0 = np.ones(x1.shape)

print(len(y))
print(len(x1))
print(len(x0))

"""
Explanation: Визуализируйте данные. По оси X — население города, по оси Y — доход грузовика.
End of explanation
"""


#w -- вектор параметров

X = np.vstack([x0, x1]).T
Y = y

w = w_init

def lossfunc(w):
    X = np.vstack([x0, x1]).T
    Y = y
    return ((1/(2*len(X))) * (X.dot(w) - Y)**2).sum()

w_init = np.array([1, 1,])

lossfunc(w_init)

def dlossfunc(w):
    X = np.vstack([x0, x1]).T
    Y = y
    return ((1/len(X)) * np.tile((X.dot(w) - Y), (2,1)).T * X).sum(axis=0)

"""
Explanation: Вспомним функцию потерь линейной регрессии:
$$
L(w) = \frac{1}{2m} \sum_{i=1}^m (h(x^i, w) - y^i)^2
$$
Здесь $h(x, w) = w^Tx = w_0 + w_1 x_1$ (предполагается, что $x_0=1$ — дополнительный признак для удобства).
$(x^i, y^i)$ — i-ый объект выборки.
Тогда правило обновления весов будет выглядеть следующим образом:
$$
w_j = w_j - lr \cdot \frac{1}{m}\sum_{i=1}^m(h(x^i, w) - y^i) x^i_j.
$$
Здесь $x^i_j$ — j-ая компонента i-ого объекта.
Определите функцию потерь и её производную. Эти функции имеют один аргумент — вектор весов $w$.
End of explanation
"""


def graddesc(w_init, lr, n): 
    #lr --- скорость обучения, шаг 
    #n --- количество итераций
    w = w_init
    l_list = []
    i_list = []
    
    for i in range (n):
        l = lossfunc(w)
        dl = dlossfunc(w)
        l_list.append(l)
        i_list.append(i)
        w = w - lr*dl
    plt.plot(i_list, l_list)
    
    return w
   


w_init = np.array([1, 1])
w = graddesc(w_init, 0.01, 10000)

print(w)

matrix([[-3.89578088],
        [ 1.19303364]])

yhat = w[0]*x0 + w[1]*x1

plt.scatter(x1, yhat)
plt.scatter(x1, y, c='red')

"""
Explanation: Напишите функцию минимизации $L(w)$ с помощью градиентного спуска, аналогичную optimize_and_plot_steps. На вход она принимает параметры обучения (темп обучения и начальное значение весов), оптимизирует итеративно функцию потерь, печатает итерации и визуализирует уменьшение функции потерь и найденное решение. Запустите функцию с постоянным темпом обучения и прокомментируйте результаты.
End of explanation
"""


# correct solution:
def softmax(x):
    """Compute softmax values for each sets of scores in x."""
    e_x = np.exp(x - np.max(x))
    return e_x / e_x.sum(axis=0)

softmax(np.array([1, 2, 3]))

"""
Explanation: Измените функцию минимизации так, чтобы темп обучения мог меняться динамически, аналогично примеру выше. Запустите функцию и прокомментируйте результаты.
Softmax
Обобщение логистической функции для многомерного случая. Функция преобразует вектор $z$ размерности $K$ в вектор $\sigma$ той же размерности, где каждая координата $\sigma_i$ полученного вектора представлена вещественным числом в интервале $[0,1]$ и сумма координат равна 1.
Координаты $\sigma_i$ вычисляются следующим образом:
${\displaystyle \sigma (z){i}={\frac {e^{z{i}}}{\displaystyle \sum {k\mathop {=} 1}^{K}e^{z{k}}}}}$

Реализуйте функцию softmax, которая на вход принимает вектор $z$, а на выходе считает от него софтмакс.
Добавьте возможность принимать на вход матрицу и считать softmax по столбцам (батч)
End of explanation
"""


