{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Векторизация\n",
    "\n",
    "## Про векторизацию\n",
    "`NumPy` оптимизирован для работы с многомерными массивами, но циклы `python` нет. В связи с этим распространен подход, называемый **векторизацией**, при котором устраняются циклы, а вместо них вызываются встроенные в `NumPy` методы. Таким образом эти циклы осуществляются внутри `C` кода. \n",
    "\n",
    "Продемонстрируем это на примере."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "N = 10 ** 7\n",
    "\n",
    "x = np.random.random(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 4.77 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "s1 = 0\n",
    "for i in range(N):\n",
    "    s1 += x[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.1 ms ± 1.24 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "s2 = x.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "`%%time` --- команда jupyter notebook измеряющая затраченное на исполнение кода в ячейке время.\n",
    "\n",
    "`%%timeit` усредняет это время по нескольким запускам.  \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " У векторизации есть три основных плюса:\n",
    " 1. Скорость: векторизованный код исполняется гораздо быстрее, чем его аналог в циклах;\n",
    " 2. Выразительность: векторизованный код больше похож на математическое выражение, что упрощает его чтение;\n",
    " 3. Количество кода: векторизованный код без циклов как правило короче, а значит в нем сложнее ошибиться."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Простые математические функции\n",
    "\n",
    "В `NumPy` есть большое количество математических функций, полный список которых можно посмотреть [здесь](https://numpy.org/doc/stable/reference/routines.math.html). Применять их можно сразу к всему массиву. Тогда выбранная функция будет применена к всем элементам массива. Все циклы выполнятся в `C` коде."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 31.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "y = np.zeros(N)\n",
    "for i in range(N):\n",
    "    y[i] = np.sin(x[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "130 ms ± 4.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "y = np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Суммы, произведения, максимум, минимум и т.п.\n",
    "\n",
    "Так же в `NumPy` реализован ряд агрегирующих методов, таких как сумма элементов массива ([ndarray.sum](https://numpy.org/doc/stable/reference/generated/numpy.sum.html#numpy.sum)), их произведение ([ndarray.prod](https://numpy.org/doc/stable/reference/generated/numpy.prod.html#numpy.prod)) и т.п. Часть из них перечислена [здесь](https://numpy.org/doc/stable/reference/routines.math.html#sums-products-differences). Методы `ndarray.min` и `ndarray.max` ссылаются на функции [np.amin](https://numpy.org/doc/stable/reference/generated/numpy.amin.html#numpy.amin) и [np.amax](https://numpy.org/doc/stable/reference/generated/numpy.amax.html#numpy.amax) соответственно.\n",
    "\n",
    "Каждая из них позволяет не только вычислять значение по всему массиву, но и вдоль обозначенной оси. Рассмотрим на примере функций `min` и `sum`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[89 85 56]\n",
      " [ 8 90 22]\n",
      " [73 44 42]]\n"
     ]
    }
   ],
   "source": [
    "x = np.random.randint(0, 100, size=(3, 3))\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n"
     ]
    }
   ],
   "source": [
    "print(x.min()) # минимум всего массива"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 8 44 22]\n"
     ]
    }
   ],
   "source": [
    "print(x.min(axis=0)) # минимум по столбцам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[56  8 42]\n"
     ]
    }
   ],
   "source": [
    "print(x.min(axis=1)) # минимум по строкам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "509\n"
     ]
    }
   ],
   "source": [
    "print(x.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[170 219 120]\n"
     ]
    }
   ],
   "source": [
    "print(x.sum(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[230 120 159]\n"
     ]
    }
   ],
   "source": [
    "print(x.sum(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Арифметика. \n",
    "\n",
    "Аналогично, всегда предпочтительнее использовать операторы `+`, `-`, `*`, `/`, `//`, `%`, `**`, `@`, `|`, `&` сразу к массивам, а не к их элементам внутри циклов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.zeros(N)\n",
    "y = np.ones(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 6.82 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "z = np.zeros(N)\n",
    "for i in range(N):\n",
    "    z[i] = x[i] + y[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46.6 ms ± 412 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "z = x + y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В ряде случаев допускается применять арифметические операторы к массивам разной формы. Строго говоря, `NumPy` может оперировать с двумя массивами, если их формы совместимы (`broadcastable`). Подробнее почитать о `broadcasting` в документации можно [здесь](https://numpy.org/doc/stable/user/theory.broadcasting.html#array-broadcasting-in-numpy) и [здесь](https://numpy.org/devdocs/user/basics.broadcasting.html).\n",
    "\n",
    "**Пример 1. Скаляр и массив.** \n",
    "\n",
    "Начнем с самого простого примера. Если в арифметическом выражении участвует массив любой формы и скаляр, то можно представить, что скаляр расширяется до массива такой же формы, у которого все элементы равны значению скаляра.\n",
    "\n",
    "![image](/_static/lecture_specific/vectorization/broadcasting_1.svg)\n",
    "\n",
    "На картинке выше массив `a` из трех элементов умножается на скаляр `b` (который можно представить в виде массива из одного элемента)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3] * [2] = [2 4 6]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "b = np.array([2])\n",
    "result = a * b\n",
    "print(f\"{a} * {b} = {result}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "В реальности никакого промежуточного массива из двоек в примере выше не создается. Это лишь удобный способ представить себе, что происходит.\n",
    "```\n",
    "\n",
    "**Общий случай.** В более общем случае, два массива совместимы, если размеры вдоль всех их осей совместимы. Размеры вдоль двух осей совместимы,    только если \n",
    "- они равны;\n",
    "- один из них равен единицы;\n",
    "При этом проверка этих размеров происходит слева направо.\n",
    "\n",
    "Рассмотрим примеры совместимых и не совместимых форм.\n",
    "\n",
    "Совместимые формы.\n",
    "\n",
    "|`A.shape`|`B.shape`|`result.shape`|Комментарий|\n",
    "| :---- | --- | --- | ---: |\n",
    "| (3,) | (1, ) | (3, ) | Размер `B` вдоль единственной оси равен 1, формы `A` и `B` совместимы |\n",
    "| (1, 3) | (3, 1) | (3, 3) | Размеры вдоль последней оси совместимы, т.к. размер `B` вдоль нее равен 1. Размеры вдоль первой оси совместимы, т.к. размер `A` вдоль неё равен 1|\n",
    "|(8, 1, 6, 1)|(7, 1, 5)|(8, 7, 5, 6)|Массив `B` можно представить в виде массива с формой (1, 7, 1, 5). Размеры вдоль всех осей совместимы, т.к. вдоль каждой оси у одного из массивов длинна равняется 1|\n",
    "|(15, 4, 3)|(15, 1, 3)|(15, 4, 3)| Размеры вдоль всех осей совместимы |\n",
    "\n",
    "Несовместимые формы.\n",
    "\n",
    "|`A.shape`|`B.shape`|Комментарий|\n",
    "| :---- | --- | ---: |\n",
    "|(3, )|(4, )| Размеры вдоль единственных осей не равны между собой, ни один из них не равен 1|\n",
    "|(2, 1)| (2, 4, 3)| Сначала сравниваются 1 и 3: всё ок. Затем сравниваются 2 и 4: несовместимы. Можно представить форму (2, 1) в виде (1, 2, 1)|\n",
    "\n",
    "**Пример 2. Матрица и вектор-строка**. \n",
    "\n",
    "Если складываются (умножаются и т.п.) двумерный и одномерный массивы, то операция возможна, если количество столбцов в матрице соответствует количеству элементов в одномерном массиве.\n",
    "\n",
    "![image](/_static/lecture_specific/vectorization/broadcasting_2.svg)\n",
    "\n",
    "В примере выше складываются массивы с формами `(4, 3)` и `(3, )`, что эквивалентно записи `(4, 3)` и `(1, 3)`. Формы совместимы и результат выходит такой, будто одномерный массив расширяется до двумерного массива с формой `(4, 3)` \"по столбцам\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  0  0]\n",
      " [10 10 10]\n",
      " [20 20 20]\n",
      " [30 30 30]] \n",
      " + \n",
      " [0 1 2] \n",
      " = \n",
      " [[ 0  1  2]\n",
      " [10 11 12]\n",
      " [20 21 22]\n",
      " [30 31 32]] \n"
     ]
    }
   ],
   "source": [
    "a = np.array([[0, 0, 0], [10, 10, 10], [20, 20, 20], [30, 30, 30]])\n",
    "b = np.array([0, 1, 2])\n",
    "result = a + b\n",
    "print(f\"{a} \\n + \\n {b} \\n = \\n {result} \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image](/_static/lecture_specific/vectorization/broadcasting_3.gif)\n",
    "\n",
    "А в этом примере формы не совместимы."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (4,3) (4,) ",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-ed51d1fba385>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0ma\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m20\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m30\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m30\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m30\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0ma\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: operands could not be broadcast together with shapes (4,3) (4,) "
     ]
    }
   ],
   "source": [
    "a = np.array([[0, 0, 0], [10, 10, 10], [20, 20, 20], [30, 30, 30]])\n",
    "b = np.array([1, 2, 3, 4])\n",
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пример 3. Матрица и вектор столбец.**\n",
    "\n",
    "Аналогично можно сложить матрицу с вектором столбцом, но для этого этот массив нужно транспонировать, т.е. привести вектор формы `(n, )` (одномерный массив из `n` элементов) к форме `(n, 1)` (матрица из `n` строк длинны 1).\n",
    "\n",
    "К сожалению, `ndarray.T` (`np.transpose`) не приведёт к желаемому результату, если массив одномерный. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4] [1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "print(b, b.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сделать это можно изменив форму или \"вставив ось\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n",
      "\n",
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n"
     ]
    }
   ],
   "source": [
    "n = b.shape[0]\n",
    "print(b.reshape((n, 1)))\n",
    "print()\n",
    "print(b[:, np.newaxis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь можно сложить матрицу со столбцом."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  1  1]\n",
      " [12 12 12]\n",
      " [23 23 23]\n",
      " [34 34 34]]\n"
     ]
    }
   ],
   "source": [
    "result = a + b[:, np.newaxis]\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Пример 3. Вектор столбец и вектор строка**\n",
    "\n",
    "Если складывать вектор столбец и вектор строку, то каждая из них как бы расширится до матрицы по строкам и по столбцам. \n",
    "\n",
    "![image](/_static/lecture_specific/vectorization/broadcasting_4.svg)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0]\n",
      " [10]\n",
      " [20]\n",
      " [30]] \n",
      " + \n",
      " [1 2 3] \n",
      "  = \n",
      " [[ 1  2  3]\n",
      " [11 12 13]\n",
      " [21 22 23]\n",
      " [31 32 33]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([0, 10, 20, 30])[:, np.newaxis]\n",
    "b = np.array([1, 2, 3])\n",
    "result = a + b\n",
    "print(f\"{a} \\n + \\n {b} \\n  = \\n {result}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}