{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explicit for-loop vs. built-ins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def index_sum_for(arr, idx):\n",
    "    sum_ = 0\n",
    "    for a, i in zip(arr, idx):\n",
    "        if i:\n",
    "            sum_ += a\n",
    "    return sum_\n",
    "\n",
    "\n",
    "def index_sum_np(arr, idx):\n",
    "    return arr[idx].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_points = int(5e5)\n",
    "array = np.random.randn((num_points))\n",
    "indices = np.random.choice([False, True], size=num_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44.7 ms ± 192 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit sum_for = index_sum_for(array, indices)\n",
    "sum_for = index_sum_for(array, indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.17 ms ± 2.25 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit sum_np = index_sum_np(array, indices)\n",
    "sum_np = index_sum_np(array, indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(sum_for, sum_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pure allocation vs. allocation + initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_points = int(1e4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.5 µs ± 33.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit out = np.empty((num_points, num_points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72.1 ms ± 626 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit out = np.full((num_points, num_points), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explicit for-loop vs. \"batched\" operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()\n",
    "A = rng.standard_normal((3, 4))\n",
    "batched_x = rng.standard_normal((1_000_000, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def apply_matrix_for(A, xs):\n",
    "    b = np.empty((xs.shape[0], A.shape[0]))\n",
    "    for i, x in enumerate(xs):\n",
    "        b[i] = A @ x\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def apply_matrix_batched(A, xs):\n",
    "    return xs @ A.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.37 s ± 20.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit b_for = apply_matrix_for(A, batched_x)\n",
    "b_for = apply_matrix_for(A, batched_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.15 ms ± 270 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit b_batched = apply_matrix_batched(A, batched_x)\n",
    "b_batched = apply_matrix_batched(A, batched_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(b_for, b_batched)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hints for Assignment 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.,  2.,  3.,  4.],\n",
       "       [ 5., 11.,  7.,  8.],\n",
       "       [ 9., 10., 16., 12.],\n",
       "       [13., 14., 15., 21.]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.set_printoptions(precision=1)\n",
    "X = np.arange(1, 17, dtype=float).reshape(4, 4) + np.eye(4) * 5\n",
    "kernel = np.array([\n",
    "    [1, 2, 1],\n",
    "    [2, 4, 2],\n",
    "    [1, 2, 1],\n",
    "]) / 16\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "def variance_naive(X):\n",
    "    M, N = X.shape\n",
    "    var_im = np.zeros_like(X)\n",
    "    for i in range(1, M - 1):\n",
    "        for j in range(1, N - 1):\n",
    "            x = X[i - 1:i + 2, j - 1:j + 2]\n",
    "            mu_x = (x * kernel).sum()\n",
    "            var_im[i, j] = (kernel * (x - mu_x) ** 2).sum()\n",
    "\n",
    "    return var_im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0. , 0. , 0. , 0. ],\n",
       "       [0. , 8.7, 9.4, 0. ],\n",
       "       [0. , 7.9, 8.7, 0. ],\n",
       "       [0. , 0. , 0. , 0. ]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "naive = variance_naive(X)\n",
    "naive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal as ss\n",
    "def variance(X):\n",
    "    # in your code this corresponds to gaussian_filter calls\n",
    "    mu_x = ss.convolve2d(X, kernel, mode='valid')\n",
    "    mu_xx = ss.convolve2d(X * X, kernel, mode='valid')\n",
    "    return np.pad(mu_xx - mu_x ** 2,  1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "smart = variance(X)\n",
    "assert np.allclose(smart, naive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}