<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Будем совсем неразумно обучаться на всем train'е, так как тогда мы переобучимся,
то есть наш алгоритм "подгониться" под закономерности, присущие только train'у,
а на реальных данных будет неистово лажать. Так что train разделим на две части
Step1:  Для каждой картинки мы хотим найти вектор $(p_0, p_1)$, вероятностей такой, что $p_i$ - вероятность того, что картинка принадлежит классу $i$ ($0$ — outdoor, $1$ — indoor).
Реализуя логистическую регрессию, мы хотим приближать вероятности к их настоящему распределению. 
Выражение выдает ответ вида  $$ W x + b, $$
где $x$ — наш вектор картинки, а результат — числовой вектор размерности $2$ с какими-то числами. Для того, чтобы эти числа стали вероятностями от $0$ до $1$, реализуем функцию 
$$
\text{softmax}(W, b, x) = \frac{e^{Wx+b}}{\sum(e^{Wx+b})},
$$
и полученные значения будут как раз давать в сумме 1, и ими мы будем приближать вероятности. 
Оценивать качество нашей модели будем с помощью кросс-энтропии, см. https
Step2:  Посчитаем среднюю квадратичную ошибку на тесте, чтобы прикинуть, что будет на Kaggle.
Step3:  Теперь применяем найденную матрицу к исследумемым данным.
<ASSISTANT_TASK:>
Python Code:
# Выделяем outdoor'ы и indoor'ы.
sample_out = sample[result[:, 0] == 1]
sample_in = sample[result[:, 1] == 1]
result_out = result[result[:, 0] == 1]
result_in = result[result[:, 1] == 1]
# Считаем размер indoor- и outdoor-частей в train'е.
train_size_in = int(sample_in.shape[0] * 0.75)
train_size_out = int(sample_out.shape[0] * 0.75)
# Разделяем  outdoor'ы и indoor'ы на обучающую и тестовую часть.
x_train_out, x_test_out = np.split(sample_out, [train_size_out])
y_train_out, y_test_out = np.split(result_out, [train_size_out])
x_train_in, x_test_in = np.split(sample_in, [train_size_in])
y_train_in, y_test_in = np.split(result_in, [train_size_in])
# Делаем общий train и test, смешивая indoor'ы и outdoor'ы.
x_train = np.vstack([x_train_in, x_train_out])
y_train = np.vstack([y_train_in, y_train_out])
x_test = np.vstack([x_test_in, x_test_out])
y_test = np.vstack([y_test_in, y_test_out])
Explanation: Будем совсем неразумно обучаться на всем train'е, так как тогда мы переобучимся,
то есть наш алгоритм "подгониться" под закономерности, присущие только train'у,
а на реальных данных будет неистово лажать. Так что train разделим на две части:
на 75% будем обучаться, а на 25% проверять, что мы лажаем не неистово.
Если возьмем первые 25% от всего train'а, то может быть несбалансированное число
outdoor'ов и indoor'ов. Поэтому для train возьмем первые 75% outdoor'ов плюс
первые 75% indoor'ов. Тогда мы сохраним пропорции outdoor:indoor таким, какое оно
во всем train'е. Будет особенно клево, если и в исследуемых данных соблюдается так же
пропорция.
End of explanation
def softmax(W, x):
    # Функция logsumexp более стабтильно вычисляет функцию экспонент, почти
    # избавляя нас от проблемы переполнения.
    p = np.dot(x, W.T)
    return np.exp(p - scm.logsumexp(p, axis=1).reshape(-1, 1))
def loss(y, softmax, W, l):
    # Формула из Википедии по ссылке выше c добавленным регуляризатором.
    return np.mean(-np.sum(y * np.log(softmax), axis=1)) + l * np.trace(W @ W.T) / (2 * y.shape[0])
# Считаем средний по всем картинкам градиент.
# Градиент у нас будет не вектор, как мы привыкли, а матрица 2x3073.
def gradients(W, x, y, l):
    p = softmax(W, x)
    grads = (p - y).T @ x + l * W
    return grads / x.shape[0]  # По максимимум матричных вычислений!
# Выбор шага по правилу Армихо из семинарского листочка.
def armijo(W, x, y, l, alpha=0.5, beta=0.5):
    s = 1
    grad = gradients(W, x, y, l)
    dW = -grad  # Направление спуска.
    loss_1 = loss(y_train, softmax(W + s * dW, x), W, l)
    loss_0 = loss(y_train, softmax(W, x), W, l)
    while loss_1 > loss_0 + alpha * s * (grad * dW).sum():
        s = beta * s
        loss_1 = loss(y_train, softmax(W + s * dW, x), W, l)
        loss_0 = loss(y_train, softmax(W, x), W, l)
    return s
def classify(x_train, x_test, y_train, y_test, iters, l):
    # Как было замечено выше, W Размера 2 на 3072, а b размера 2, но мы приписываем b к W.
    W = np.zeros((2, 3072))
    b = np.zeros(2)
    # Для приписывания запишем b как вектор столбец и воспользуемся функцией hstack.
    b = b.reshape(b.size, 1)
    W = np.hstack([W, b])
    # Соответственно, нужно поменять x_train и x_test, добавив по 1 снизу.
    fictious = np.ones((x_train.shape[0], 1))
    x_train = np.hstack([x_train, fictious])
    fictious = np.ones((x_test.shape[0], 1))
    x_test = np.hstack([x_test, fictious])
    # Будем записывать потери на каждом шаге спуска.
    losses_train = [loss(y_train, softmax(W, x_train), W, l)]
    losses_test = [loss(y_test, softmax(W, x_test), W, l)]
    # Собственно, сам спуск.
    for i in tqdm.tqdm(np.arange(iters)):
        # Именно так - в Армихо подставляется alpha = l, а l = 0!
        # Потому что я накосячил и не заметил! =)
        eta = armijo(W, x_train, y_train, 0, l)
        W = W - eta * gradients(W, x_train, y_train, l)
        losses_train.append(loss(y_train, softmax(W, x_train), W, l))
        losses_test.append(loss(y_test, softmax(W, x_test), W, l))
    # На выходе имеется оптимальное значение W и массивы потерь.
    return W, losses_train, losses_test
l = 0.04  # Сработает лучше, чем вообще без регуляризатора (l = 0).
# Нам хватит и 100 итераций, переобучение начинается достаточно быстро.
W, losses_train, losses_test = classify(x_train, x_test, y_train, y_test, 100, l)
plt.plot(losses_train, color='green', label='train')
plt.plot(losses_test, color='red', label='test')
plt.xlabel('Gradient descent iteration')
plt.ylabel('Loss')
plt.legend()
plt.show()
iters = np.argmin(losses_test)  # На этой итиреации ошибка на тесте минимальна.
# Делаем столько итераций.
W, losses_train, losses_test = classify(x_train, x_test, y_train, y_test, iters, l)
Explanation: Для каждой картинки мы хотим найти вектор $(p_0, p_1)$, вероятностей такой, что $p_i$ - вероятность того, что картинка принадлежит классу $i$ ($0$ — outdoor, $1$ — indoor).
Реализуя логистическую регрессию, мы хотим приближать вероятности к их настоящему распределению. 
Выражение выдает ответ вида  $$ W x + b, $$
где $x$ — наш вектор картинки, а результат — числовой вектор размерности $2$ с какими-то числами. Для того, чтобы эти числа стали вероятностями от $0$ до $1$, реализуем функцию 
$$
\text{softmax}(W, b, x) = \frac{e^{Wx+b}}{\sum(e^{Wx+b})},
$$
и полученные значения будут как раз давать в сумме 1, и ими мы будем приближать вероятности. 
Оценивать качество нашей модели будем с помощью кросс-энтропии, см. https://en.wikipedia.org/wiki/Cross_entropy.
Сначала поймем, что $x$ - вектор размерности 3072, $W$ - матрица 2 на 3072, $b$ - вектор размерности 2.
Положим $x'i = x_i$ для $ i \leqslant 3072 $ и $x'{3073} = 1$. Получили вектор $x'$ размерности 3073. Положим $W'{i,j} = W{i,j}$ для $ i \leqslant 2, j \leqslant 3073$ и $W'_{i,3073}=b_i$ для $ i \leqslant 2 $.
Таким образом, к вектору $x$ просто дописали 1, а к матрице $W$ просто приписали вектор $b$ справа.
Заметим теперь, что в точности верно равенство: $Wx+b=W'x'$. Теперь забьем на вектор $b$ и будем считать, что у нас есть матрица 10 на 3073, элементы которой надо оценить. Далее везде считаем $W' = W$ и $x' = x$.
Градиентный спуск считается по формуле: $W_{k+1} = W_k - \eta_k \nabla L(W_k)$, где $\eta_k$ — шаг, а $L$ — функция $\text{loss}$. Значит, нам надо посчитать градиент функции $L$, то есь найти ее частные производные по всем 6146 переменным.
Вспомним, как определяется $L$. Обозначим через $y$ вектор вида $(1, 0)$ либо $(0, 1)$, где 1 на $k$-м месте, где $k - 1$ — тип исследуемой картинки. Размерность $y$ равна 2. Сам вектор $y$ олицетворяет ответ для данной картинки.
Тогда
$$ L(W) = -y_1 \ln \frac{e^{(Wx)1}}{e^{(Wx)_1} + e^{(Wx)_2}} -y{2} \ln \frac{e^{(Wx){2}}}{e^{(Wx)_1} + e^{(Wx)_2}} + \frac{\lambda}{2} \sum{i=1}^{2} \sum_{j=1}^{3073} W_{i,j}^2. $$
Последняя сумма — так называемый регуляризатор. Если у нас много признаков (у нас их 6146), то при логистической регресии может возникнуть переобучение. Добавляя все параметры в $\text{loss}$, мы не сможем получить неестественного результата, когда какие-то параметры очень маленькие, а какие-то очень большие, потому что большие будут сильно увеличивать регуляризатор, а функция минимизируется. Таким образом, более вероятно получение подходящего результата.
Это описано в курсе Machine Learning by Stanford University во втором уроке третьей недели. Ссылка: https://www.coursera.org/learn/machine-learning/lecture/4BHEy/regularized-logistic-regression.
Теперь найдем производную по $W_{i,j}$: $$
\frac{dL(W)}{dW_{i,j}} =
-y_1 \frac{e^{(Wx)1} + e^{(Wx)_2}}{e^{(Wx)_1}} \cdot
\frac{-e^{(Wx)_1} e^{(Wx)_i} x_j}
{e^{(Wx)_1} + e^{(Wx)_2}}
-y{2} \frac{e^{(Wx)1} + e^{(Wx)_2}}{e^{(Wx){2}}} \cdot
\frac{-e^{(Wx){2}} e^{(Wx)_i} x_j}
{e^{(Wx)_1} + e^{(Wx)_2}} -\
- y_i \frac{e^{(Wx)_1} + e^{(Wx)_2}}{e^{(Wx)_i}} \cdot
\frac{e^{(Wx)_i} x_j (e^{(Wx)_1} + e^{(Wx)_2})}
{(e^{(Wx)_1} + e^{(Wx)_2})^2}
+ \lambda W{i,j}. $$
Упростим немного: $$
\frac{dL(W)}{dW_{i,j}} =
\frac{ x_j e^{(Wx)i} (y_1 + y_2) }
{e^{(Wx)_1} + e^{(Wx)_2}}
-y_i x_j
+ \lambda W{i,j}. $$
Упрощая еще сильнее, приходим к окончательному ответу:  $$
\frac{dL(W)}{dW_{i,j}} =\left( \frac{e^{(Wx)i}}{e^{(Wx)_1} + e^{(Wx)_2}} - y_i \right) x_j
+ \lambda W{i,j}.
$$
Соответственно, если $j = 3073$, то есть дифференцируем по переменным $ W_{1, 3073} = b_1, \ldots, W_{2, 3073} = b_2$, то коэффициент перед скобкой просто 1.
Перейдем к реализации.
End of explanation
# Добавляем 1 к выборке.
nx_test = np.hstack([x_test, np.ones(x_test.shape[0]).reshape(x_test.shape[0], 1)])
probabilities = softmax(W, nx_test)  # Считаем вероятности.
recognized = np.argmax(probabilities, axis=1)  # Что распознано.
answers = np.argmax(y_test, axis=1)  # Правильные ответы.
np.sqrt(np.mean((recognized - answers) ** 2))  # Собственно, ошибка.
Explanation: Посчитаем среднюю квадратичную ошибку на тесте, чтобы прикинуть, что будет на Kaggle.
End of explanation
# Добавляем 1 к выборке.
ntest = np.hstack([test, np.ones(test.shape[0]).reshape(test.shape[0], 1)])
probabilities = softmax(W, ntest)  # Считаем вероятности.
ress = np.argmax(probabilities, axis=1).reshape(-1, 1)  # Что распознано.
# Осталось загнать все в табличку, чтобы ее записать в csv.
ids = np.arange(ress.size).reshape(-1, 1)
submit = np.hstack([ids, ress])
# Заполняем csv-шник.
import csv
with open('submission.csv', 'w', newline='') as csvfile:
    submission = csv.writer(csvfile, delimiter=',')
    submission.writerow(['id', 'res'])
    submission.writerows(submit)
Explanation: Теперь применяем найденную матрицу к исследумемым данным.
End of explanation
<END_TASK>