{ "cells": [ { "cell_type": "markdown", "id": "408032c3", "metadata": {}, "source": [ "# STAT 540 — Programming Basics in Python\n", "## Building a Grading Machine\n", "\n", "**Topics:** `if` / `elif` / `else` · vectorized conditionals · `for` loops · functions · `map` / comprehensions / `apply`\n", "\n", "This notebook is the Python companion to the R lecture `programming-basics_9_09_2026.R`.\n", "Every R idea has a Python counterpart; where they differ, a **R vs. Python** note points it out." ] }, { "cell_type": "code", "execution_count": 3, "id": "50fb448e", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from scipy import stats\n", "\n", "np.set_printoptions(precision=4, suppress=True)" ] }, { "cell_type": "markdown", "id": "37c3406f", "metadata": {}, "source": [ "---\n", "# 1. `if` / `else`" ] }, { "cell_type": "markdown", "id": "156b8a11", "metadata": {}, "source": [ "### The reciprocal example\n", "\n", "```r\n", "# R\n", "a <- 0\n", "if (a != 0) {\n", " print(1/a)\n", "} else {\n", " print(\"No reciprocal for 0.\")\n", "}\n", "```\n", "\n", "Python uses **indentation** instead of braces, and a **colon** after the condition." ] }, { "cell_type": "code", "execution_count": 1, "id": "da70e435", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No reciprocal for 0.\n" ] } ], "source": [ "a = 0\n", "\n", "if a != 0:\n", " print(1 / a)\n", "else:\n", " print(\"No reciprocal for 0.\")" ] }, { "cell_type": "markdown", "id": "59cfaf73", "metadata": {}, "source": [ "> **R vs. Python**\n", "> | R | Python |\n", "> |---|---|\n", "> | `if (cond) { ... } else { ... }` | `if cond:` … `else:` |\n", "> | braces `{ }` delimit blocks | **indentation** (4 spaces) delimits blocks |\n", "> | `!=`, `==`, `>=` | same |\n", "> | &&, ||, ! | `and`, `or`, `not` |" ] }, { "cell_type": "markdown", "id": "aa5f7873", "metadata": {}, "source": [ "### The murders data\n", "\n", "`dslabs::murders` is an R package dataset. In Python we read the same data with **pandas**." ] }, { "cell_type": "code", "execution_count": 7, "id": "cf109cb1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stateabbregionpopulationtotal
0AlabamaALSouth4779736135
1AlaskaAKWest71023119
2ArizonaAZWest6392017232
3ArkansasARSouth291591893
4CaliforniaCAWest372539561257
\n", "
" ], "text/plain": [ " state abb region population total\n", "0 Alabama AL South 4779736 135\n", "1 Alaska AK West 710231 19\n", "2 Arizona AZ West 6392017 232\n", "3 Arkansas AR South 2915918 93\n", "4 California CA West 37253956 1257" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Reads the local copy if it is sitting next to this notebook,\n", "# otherwise pulls it straight from the dslabs GitHub repo.\n", "try:\n", " murders = pd.read_csv(\"murders.csv\")\n", "except FileNotFoundError:\n", " murders = pd.read_csv(\n", " \"https://raw.githubusercontent.com/rafalab/dslabs/master/inst/extdata/murders.csv\"\n", " )\n", "\n", "murders.head()" ] }, { "cell_type": "code", "execution_count": 13, "id": "10767a09-0464-4947-b629-340b5a0ba5ec", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stateabbregionpopulationtotal
0AlabamaALSouth4779736135
1AlaskaAKWest71023119
2ArizonaAZWest6392017232
3ArkansasARSouth291591893
4CaliforniaCAWest372539561257
\n", "
" ], "text/plain": [ " state abb region population total\n", "0 Alabama AL South 4779736 135\n", "1 Alaska AK West 710231 19\n", "2 Arizona AZ West 6392017 232\n", "3 Arkansas AR South 2915918 93\n", "4 California CA West 37253956 1257" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "murders.iloc[:5]" ] }, { "cell_type": "code", "execution_count": 15, "id": "19b4f034-0fac-4ef7-83fe-30b2e6bbb435", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stateabbregion
0AlabamaALSouth
1AlaskaAKWest
2ArizonaAZWest
3ArkansasARSouth
4CaliforniaCAWest
\n", "
" ], "text/plain": [ " state abb region\n", "0 Alabama AL South\n", "1 Alaska AK West\n", "2 Arizona AZ West\n", "3 Arkansas AR South\n", "4 California CA West" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "murders.iloc[:5, :3]" ] }, { "cell_type": "code", "execution_count": 10, "id": "6316a111", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 2.824424\n", "1 2.675186\n", "2 3.629527\n", "3 3.189390\n", "4 3.374138\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# R: murder_rate <- murders$total / murders$population * 100000\n", "murder_rate = murders[\"total\"] / murders[\"population\"] * 100000\n", "murder_rate.head()" ] }, { "cell_type": "markdown", "id": "f9f67553", "metadata": {}, "source": [ "### `which.min()` → `argmin()`\n", "\n", "> **R vs. Python** R indexes from **1**; Python indexes from **0**.\n", "> `which.min(x)` returns a 1-based position, `np.argmin(x)` / `Series.argmin()` returns a 0-based one." ] }, { "cell_type": "code", "execution_count": 16, "id": "ad34ce1b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "position: 45\n", "state: Vermont\n", "rate: 0.31962105727449536\n" ] } ], "source": [ "ind = murder_rate.argmin() # 0-based position of the smallest rate\n", "print(\"position:\", ind)\n", "print(\"state:\", murders[\"state\"][ind])\n", "print(\"rate:\", murder_rate[ind])" ] }, { "cell_type": "code", "execution_count": 17, "id": "63e8b092", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Vermont\n" ] } ], "source": [ "if murder_rate[ind] < 0.5:\n", " print(murders[\"state\"][ind])\n", "else:\n", " print(\"No state has murder rate that low\")" ] }, { "cell_type": "code", "execution_count": 18, "id": "a30e21a7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No state has a murder rate that low.\n" ] } ], "source": [ "if murder_rate[ind] < 0.25:\n", " print(murders[\"state\"][ind])\n", "else:\n", " print(\"No state has a murder rate that low.\")" ] }, { "cell_type": "markdown", "id": "d4ff4c13", "metadata": {}, "source": [ "---\n", "# 2. Exercise — Build a Grading Machine (Pass / Fail)" ] }, { "cell_type": "code", "execution_count": 19, "id": "93bbcaf5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[55, 82, 91, 40, 67, 78, 100, 62, 35, 88]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scores = [55, 82, 91, 40, 67, 78, 100, 62, 35, 88]\n", "scores" ] }, { "cell_type": "markdown", "id": "d912a2e1", "metadata": {}, "source": [ "**Task.** Write an `if`/`else` statement that reports **PASS** (score ≥ 60) or **FAIL** for the *second* score.\n", "\n", "> **Careful!** The R version used `scores[2]`. In Python the second element is `scores[1]`." ] }, { "cell_type": "code", "execution_count": 20, "id": "258f3d29", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PASS\n" ] } ], "source": [ "if scores[1] >= 60: # scores[1] is 82 -- the 2nd student\n", " print(\"PASS\")\n", "else:\n", " print(\"FAIL\")" ] }, { "cell_type": "markdown", "id": "5c9f3b70", "metadata": {}, "source": [ "---\n", "# 3. Exercise 2 — Letter grades with `elif`\n", "\n", "```\n", "score >= 90 -> A\n", "score >= 80 -> B\n", "score >= 70 -> C\n", "otherwise -> F\n", "```\n", "\n", "> **R vs. Python** R's `else if` becomes a single keyword in Python: **`elif`**." ] }, { "cell_type": "code", "execution_count": 21, "id": "d1d047df", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "F\n" ] } ], "source": [ "if scores[8] >= 90: # scores[8] is 35 -- the 9th student\n", " print(\"A\")\n", "elif scores[8] >= 80:\n", " print(\"B\")\n", "elif scores[8] >= 70:\n", " print(\"C\")\n", "else:\n", " print(\"F\")" ] }, { "cell_type": "markdown", "id": "8eac9a29", "metadata": {}, "source": [ "---\n", "# 4. Vectorized conditionals: `ifelse()` → `np.where()`\n", "\n", "R's `ifelse()` does two different jobs. Python splits them:\n", "\n", "| R | Python (one value) | Python (whole vector) |\n", "|---|---|---|\n", "| `ifelse(a > 0, 1/a, NA)` | `1/a if a > 0 else np.nan` | `np.where(a > 0, 1/a, np.nan)` |" ] }, { "cell_type": "code", "execution_count": 22, "id": "0b8d1ba1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "nan" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# --- single value: the conditional expression ---\n", "a = 0\n", "result = 1 / a if a > 0 else np.nan\n", "result" ] }, { "cell_type": "code", "execution_count": 23, "id": "0471493f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([nan, 1. , 0.5, nan, 0.2])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# --- whole array: np.where ---\n", "a = np.array([0, 1, 2, -4, 5], dtype=float)\n", "\n", "with np.errstate(divide=\"ignore\"): # silence the 1/0 warning\n", " result = np.where(a > 0, 1 / a, np.nan)\n", "\n", "result" ] }, { "cell_type": "code", "execution_count": 24, "id": "15ee3780", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['MALE', 'FEMALE', 'FEMALE', 'MALE'], dtype=' **R vs. Python** `1:5` in R includes 5. Python's `range(1, 6)` **stops before** 6." ] }, { "cell_type": "code", "execution_count": 27, "id": "a12182a9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "2\n", "3\n", "4\n", "5\n" ] } ], "source": [ "for i in range(1, 6):\n", " print(i)" ] }, { "cell_type": "markdown", "id": "c9b1c24f", "metadata": {}, "source": [ "### Building up a vector inside a loop\n", "\n", "R:\n", "```r\n", "m <- 25\n", "s_n <- vector(length = m)\n", "for (n in 1:m) {\n", " s_n[n] <- sum(1:n)\n", "}\n", "```" ] }, { "cell_type": "code", "execution_count": 28, "id": "1d9a7e94", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1., 3., 6., 10., 15., 21., 28., 36., 45., 55., 66.,\n", " 78., 91., 105., 120., 136., 153., 171., 190., 210., 231., 253.,\n", " 276., 300., 325.])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = 25\n", "s_n = np.zeros(m) # an empty vector of the right size\n", "\n", "for n in range(1, m + 1):\n", " s_n[n - 1] = sum(range(1, n + 1)) # careful: position n-1 holds the n-th answer\n", "\n", "s_n" ] }, { "cell_type": "code", "execution_count": null, "id": "a3147ee0", "metadata": {}, "outputs": [], "source": [ "# Avoid the n-1 bookkeeping with enumerate()\n", "s_n = np.zeros(m)\n", "for i, n in enumerate(range(1, m + 1)):\n", " s_n[i] = np.arange(1, n + 1).sum()\n", "\n", "s_n[:5]" ] }, { "cell_type": "code", "execution_count": 31, "id": "9da24623-e27d-424b-b726-8438aff26631", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 1\n", "1 2\n", "2 3\n", "3 4\n", "4 5\n", "5 6\n", "6 7\n", "7 8\n", "8 9\n", "9 10\n", "10 11\n", "11 12\n", "12 13\n", "13 14\n", "14 15\n", "15 16\n", "16 17\n", "17 18\n", "18 19\n", "19 20\n", "20 21\n", "21 22\n", "22 23\n", "23 24\n", "24 25\n" ] } ], "source": [ "m=25\n", "for i, n in enumerate(range(1, m + 1)):\n", " print(i, n)" ] }, { "cell_type": "markdown", "id": "666bbe82", "metadata": {}, "source": [ "---\n", "# 6. Exercise — a `for` loop over all the scores\n", "\n", "Return the letter grade for **every** student." ] }, { "cell_type": "code", "execution_count": 32, "id": "6557a68a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "F\n", "B\n", "A\n", "F\n", "F\n", "C\n", "A\n", "F\n", "F\n", "B\n" ] } ], "source": [ "scores = [55, 82, 91, 40, 67, 78, 100, 62, 35, 88]\n", "\n", "for i in scores:\n", " if i >= 90:\n", " print(\"A\")\n", " elif i >= 80:\n", " print(\"B\")\n", " elif i >= 70:\n", " print(\"C\")\n", " else:\n", " print(\"F\")" ] }, { "cell_type": "markdown", "id": "b6731080", "metadata": {}, "source": [ "---\n", "# 7. Functions\n", "\n", "```r\n", "# R\n", "newfunction <- function(x, y){\n", " ans <- x^2 + y^2\n", " return(ans)\n", "}\n", "```\n", "\n", "> **R vs. Python**\n", "> `function(x, y)` → `def newfunction(x, y):` · `^` → `**` · `return(ans)` → `return ans`" ] }, { "cell_type": "code", "execution_count": 33, "id": "3e53616b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def newfunction(x, y):\n", " ans = x**2 + y**2\n", " return ans\n", "\n", "\n", "newfunction(3, 4)" ] }, { "cell_type": "code", "execution_count": 34, "id": "e1d5994c", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE1CAYAAAAoBaPAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOxpJREFUeJzt3QlYVGXbB/A/+74IyKYgCCgqbqm4ay6Z2mJWbrnnW7lVWmr6fqX5ZmllZZmZWbmbqbmUa+a+4S7ihsomogiI7LLzXc8zQFpaiMOcMzP/X9dcnJlG5uEcmHue7b5NSkpKSkBERESqZKp0A4iIiOjBGKiJiIhUjIGaiIhIxRioiYiIVIyBmoiISMUYqImIiFSMgZqIiEjFzKHHiouLcf36dTg4OMDExETp5hAREVWISGGSmZkJb29vmJqaGm6gFkHax8dH6WYQERFVSnx8PGrWrGm4gVr0pMt+UEdHR6WbQ0REVCEZGRmyo1kWxww2UJcNd4sgzUBNRET6piLTtlxMRkREpGIM1ERERCrGQE1ERKRiDNREREQqxkBNRESkYgzUREREKsZA/Re5BUXKXAkiIlKtvMIiFBWXKPLaDNSlIhMzMeC7MLy8+JgiF4KIiNRr6aE4dP5sD34Nv67z19brhCfaZGdlhqOxqfIT07nr6Wjg7aR0k4iISAUKi4qx+FAsEtLu4E5+oc5fnz3qUjWr2aJHiKc8/uFAjM4vBBERqdPv52/KIO1iZ4leTWro/PUZqO/yn/a15dffwq/jZkauzi8GERGpz4+lnbdBLX1hbWGm89dnoL5LEx9ntPCrhoKiEiw9HKvzi0FEROpyOj4Nx+Nuw8LMBINa1VKkDQzUfzGinaZXveLIVeQoMBdBRETqsXB/tPz6bOMacHe0VqQNDNR/8UR9D/i62CItpwC/nExQ5KIQEZHy4lNzsDXihjz+T3t/xdrBQP0XZqYmeLmtX/m8RLFC++aIiEhZYqW3CAHtg9xQz8tRsXYwUN9Hn+Y+cLA2R0xKNnZdTNL9VSEiIkVl5Bbg52Px8nhEO+V60wID9X3YWZnjpZa+8vj7A5r5CSIiMh6rjl5FVl4hgtzt0bFOdUXbwkD9AMPa+MHc1ARh0ak4m5Cu26tCRESKKRAJTg7Gls9Nm5iYKHo1GKgfwMvJBk818pLHTIBCRGQ8tkTcwPX0XLjZK5Pg5K8YqP/Bf9r9mQAlMZ0JUIiIDF1JSQm+369JcDKktZ8iCU7+ioH6HzSs6YSW/i4oLC7BEiZAISIyeEdiUhGRkA4rc1MMLF2rpDQG6gqmFV0RFofsPCZAISIyZN+X9qZfaFYTrvZWUAMG6n/RJdgdfq62yMgtxC8nr+nmqhARkc5FJ2dh58WbqtiSdTcG6n87QaYm5RdMLCpTqnA4ERFVLfEeX1Ki6aAFVLeHWjBQV4AYAnGysUDcrRzsOJ9Y9VeFiIh0KjU7v3zUtGzKUy0YqCvA1tIcQ1prqqbM3xstVwUSEZHhWB4Wh9yCYjTwdkSr2i5QEwbqChraxk+uAgyPT5OrAomIyDDkFhTJvN7Cqx1qK57g5K8YqCvIzd4KfZrXlMcL9kZV5TUhIiIdWnM8Xg5916xmg6caahJdqQkD9UMmQDE1AXZHJuNiYkbVXRUiItKJwqJifFdac/qV9rVhbqa+sKi+FqmYn5sdeoRoPm19t4/FOoiI9N3Ws4mIT72DarYW5aOmaqNooC4qKsJ7770Hf39/2NjYICAgAB988IGqF2u91lGzGvDX09eRkHZH6eYQEVEliVjzbelUpliHJBYOq5Gigfrjjz/G/Pnz8fXXX+PChQvy/ieffIK5c+dCrRrVdEabAFeZVvTHA5oMNkREpH8OXEnBuesZsLEww9DWflArRQP1oUOH0KtXLzz11FPw8/PDiy++iG7duuHo0aP3fX5eXh4yMjLuuSnhtY4B8utPR68iPadAkTYQEdGjKetN92vhg2p2llArRQN1mzZtsHPnTly6dEneDw8Px4EDB9CjR4/7Pn/mzJlwcnIqv/n4+EAJHYLcUM/LETn5RVgWplnST0RE+iPiWjoOXrkFs7uyT6qVooF68uTJ6N+/P4KDg2FhYYGmTZti3LhxGDhw4H2fP2XKFKSnp5ff4uPjoQSxx25k6Vy12Hsn9uAREZH++Hafpjf9TCMv+LjYQs0UDdSrV6/GihUrsHLlSpw8eRJLlizB7Nmz5df7sbKygqOj4z03pfRs6IUazjZIyfoz7RwREalf3K1sbI24cc9UppopGqgnTpxY3qtu2LAhBg8ejPHjx8shbrWzMDPFf9prhksW7otmsQ4iIj2xcH80RH2ljnWqy2lMtVM0UOfk5MDU9N4mmJmZobi4GPpALEBwtrVA7K0cbD/HYh1ERGqXkpWHNcc1o6Aj9aA3rXigfuaZZ/Dhhx9i8+bNiI2Nxfr16/H555+jd+/e0J9iHZol/fP3RKl6/zcREQGLD8Yir7AYjX2cVVd8Q5WBWuyXFluyRo8ejXr16mHChAl47bXXZNITfTGsjZ/cgxeRkI59l1OUbg4RET1AVl4hlh7W7NQZ1VF9xTdUGagdHBwwZ84cxMXF4c6dO4iKisKMGTNgaane/Wx/5WJniZda+srjebuvKN0cIiL6h1KWGbmFqO1mhyfqe0JfMNe3FohE7pZmpjgak4pjsSyBSUSkNrkFRfh+vyab5KjHA+T+aX3BQK0Fnk7WeKGZJpn717vYqyYiUpufj8XLhWRiW+1zTWtAnzBQa8mojgGyBObeS8k4m5CurW9LRESPKL+wGAtK04WKZFVie60+0a/Wqpivqy2ebewtjzlXTUSkHhtOJeB6ei6qO1ihT3NlUk8/CgZqLRrdKVB+3XYuEVeSMrX5rYmIqBKKikswv7Q3/Up7f1hbmEHfMFBrUR0PBzzZwANiO/U3ezS/GEREpJzNETcQk5Itk1MNbFlLLy8FA7WWjSntVW88fR3xqTna/vZERFRBxcUl+KZ02+zwNv6wszKHPmKg1rJGNZ3RPshNDreU1TolIiLd23kxCRcTM2FvZS6TU+krBuoqMLa0Vy3yySZl5FbFSxAR0T8QKZ2/Lu1ND2pVC062FtBXDNRVINTfBc1rVUN+UbGs0kJERLp18MothMenwcrcFCPaaSod6isG6iog8seO6azpVa84chW3s/Or4mWIiOgBvt59WX4dEOort2XpMwbqKvJ4nepo4O2InPwifH+AvWoiIl05EZeKsOhUWJiZ4NUOtfX+xDNQV2Gv+o0uQfJ4yaE4pOWwV01EpAtf7dTMTT/ftCa8nW30/qQzUFehbvU9UM/LUZZW++GAJhk8ERFVnVNXb8tUzqLoxuhOAQZxqhmoq7hX/WYXzVz1ooOx7FUTEVWxL3dq5qafb1oDtVztDOJ8M1BXsW71PRHs6SB71T+yV01EVKW96T2Rmt702NIFvYaAgbqqT7Cp6FUHlfeq03MKqvoliYiM0pcG2JsWGKh14MkGnqjr4YBM0as+yLlqIiJtOx2fZpC9aYGBWhcn2fTPFeAiUKffYa+aiEibvvzjkvza28B60wIDtY70CPFEHQ97ZOYWYhF71UREWu1N7y7rTZemcDYkDNRK9KoPsFdNRFQVvWk/N8PqTQsM1DrUM8QLQe72yMgtxOKDsbp8aSIigxRu4L1pgYFalyfb1ASvl/aqfzgQjYxczlUTEWljpfdzTQyzNy0wUOvYUw29EFDdTvaql7BXTUT0SL3pXReTZG/6dQNb6X03BmodM7trrlqUwOQKcCKiyvnSCHrTAgO1Ap5u5F0+V/0D61UTEVVqpfcuI+hNCwzUChC/WG93qyOPRbGOW1l5SjSDiEhvffZ7pEGv9L4bA7WC2cpEvers/CIs2Md61UREFRUWfQv7L6fIetNlKZoNGQO1gpW1JnSrK4+XHIpFUkauUk0hItIbJSUlmL1d05vu38IXPi62MHQM1Ap6vG51NKtVDXmFxZi3W1PonIiIHmzPpWQcj7sNK3NTg8vp/SAM1Ar3qsvmqlcevYprt3OUbA4Rkep705+Vzk0PaV0LHo7WMAYM1AprE+CGNgGuKCgqwdyd7FUTET3ItrOJOJuQATtLM4x63Dh60wIDtQq8XTpXvfbkNcSkZCvdHCIi1SkqLsFnOzQ5vUe084eLnSWMBQO1Coh56s7B7vIXsSy5PBER/enX8ARcScqCo7U5RrSvbVSnhoFaJd56QjNXvTH8Oi7dzFS6OUREqlFQVIwvdmiykL3WMQBONhYwJpUO1FevXsX+/fuxfft2nDx5Enl5TNrxKEJqOKFnQ0+UlABflA7vEBERsOb4NVxNzYGbvSWGt/UzulPyUIE6NjYW77zzDmrVqgV/f3907NgRPXr0QPPmzeHk5IQnnngCa9asQXFxcdW12ICN71oHJibA1rOJiLiWrnRziIgUl1tQhLm7NL3p0Y8HwtbSHMamwoH6jTfeQOPGjRETE4MZM2bg/PnzSE9PR35+PhITE7Flyxa0a9cOU6dORaNGjXDs2LGqbbkBCvJwQO8mNeTxJ9svKt0cIiLFrThyFTfSc+HlZI2XWvrCGFU4UNvZ2SE6OhqrV6/G4MGDUbduXTg4OMDc3Bzu7u7o3Lkzpk2bhgsXLmD27NmIj4+v0PdNSEjAoEGD4OrqChsbGzRs2BDHjx+HsRr/RB1YmpnK9Hj7Lycr3RwiIsVk5haUJ4N6vXMQrC3MjPJqVHgMYebMmRX+pt27d6/Q827fvo22bduiU6dO2Lp1K6pXr47Lly+jWrVqMFYiHd6gVrXw48EYfLztItoGuMHU1ETpZhER6dx3+6KRmp2P2m526NO8ptFegUotJrt48cHDsmJxWUV9/PHH8PHxwaJFixAaGirnvbt164aAgAAYM5EWz97KXG7s3xRxQ+nmEBHpnKh/8P3+GHk8qXtdWJgZ7yalSv3kjz32GObNm3fPY2LV99ixY9GrV68Kf59ff/1VLkTr06ePHD5v2rQpFi5c+MDni9fIyMi452aIxEb+1zpo9gmKdHn5hVycR0TG5cudl3GnoAhNfJxltUFjVqlAvXjxYrlorGfPnrh58yZOnz4tg+wff/wht2xVlJjznj9/PoKCgmRPfNSoUXLR2pIlSx44/C5Wl5fdRG/cUI1o7w83eyvE3crBqmNXlW4OEZHORCdnYdUxzTqnKT2CZV0EY2ZSIrKcV8K1a9cwfPhwnDp1CtnZ2Rg2bBg+++wz2NpWvOSYpaWl7FEfOnSo/DERqMWK8cOHD9+3R333fm3RoxbBWqw+d3R0hKFZFhaH9zaclXsH90zsJIfDiYgM3egVJ7AlIhFdgt3xw7AWMEQifokOZ0Xi1yMN+outWUVFRfLm5eUFa+uHq2Qi/k39+vXveaxevXoymcr9WFlZyR/o7psh69/CB36utkjJysf3+6OVbg4RUZU7dfW2DNKiEz2pezDPeGUD9apVq+Q2KvFp4NKlS9i8eTO+++47tG/fXg5nV5RY8R0ZqSlZVkZ8P5FQhSAXT0x4UlOwY+G+aKRkMfsbERkuMcA7c6tmsfILj9VEXU8HpZukv4F6xIgR+Oijj+RiMLGlSmQki4iIQI0aNdCkSZMKf5/x48cjLCxMfq8rV65g5cqVMuCPGTOmMs0ySD1DvNCophOy84vw9S6WwSQiw7UnMhlHY1JhaW5aXv+AKhmoRW5vsfDrbmLvs0iG8tfV4P+kRYsWWL9+PX766SeEhITggw8+wJw5czBw4EBem1JiD/Xk0uGfFUfiEHeLZTCJyPCI6oEid4QwvI0fvJ1tlG6S/i8m07fJeH03+IcjMlvZM429MXdAU6WbQ0SkVb+cuIa314TLMpb7JnWCs61h15vOqIrFZLNmzcKdO3cq9NwjR47IeWvSnne6iy0KwG/h13E6Po2nlogMqvDG56VVA0d3CjT4IP2wKhyoRREOX19fjB49Wqb7TE7+Mw91YWEhzpw5g2+++QZt2rRBv379ZB5w0m4ZzN5NNQU7Zmw6LxddEBEZgkUHY5GQdkcW3hjWxvjKWGotUC9dulQmNCkoKMBLL70ET09PuQ9aBGSxbUokPPnxxx8xZMgQmWK0Q4cOFf3WVEETn6wLawtTHI+7jW1nE3neiEjvid0s35QW3tC8xxln4Q2tz1GLetOiBx0XFyeHw93c3ORqb/FVl4xpjrqMGB76audl+LrYYsdbHWBlzl9qItJf726IwPKwq2hYwwkbx7Q1miJEGQ8RvyqV6srU1FQG5ofZikXaIXKA/3T0Kq6m5mDpoTi8UpoTnIhI31y6mYmVRzQJrt59qp7RBOmHVenMZFFRUXj33XcxYMAAJCUlycfE3PW5c+cq+y2pAuyszDGxmyYJyle7LssScERE+uijLRdQXAI82cADLWu7Kt0cwwrUe/fulZnJxOrudevWISsrSz4eHh6OadOmabuN9BcvNKuJYE8HZOYWymFwIiJ9s/dSskxwYmFmgsk96indHMML1JMnT8aMGTOwY8cOuaCsTOfOnWWmMapaZqYmePcpTY705WFxiErWfFAiItIHhUXF+HDzeXk8pLUf/N3slG6S4QVqkS60d+/ef3tc1JROSUnRRrvoX7QLcpOVZQqLSzBziyabDxGRPlh9/Bou3cyCs60F3ugcpHRzDDNQOzs748aNG397XJS8FPm+STem9Kwne9d/XLiJQ1H8gERE6peVV4jPd2iKMb3ZJQhOthZKN8kwA3X//v3xzjvvIDFRlCIzkdu1Dh48iAkTJsh91KQbge72GNjSVx5/uPkCisWqDCIiFZu/54os3SuGuwe2ZKXEKgvUotpVcHAwfHx85EIyUVNaJDgRWcnESnDSHfGJ1MHKHOeuZ2DtyWs89USkWvGpOfh+f4w8ntIjWFbJon9XqbMkFpAtXLhQbtHatGkTli9fLrORLVu2DGZmTMChS672Vni9S6A8/mRbJDJyC3T6+kRED7MdK6+wGK1ru+KJ+h48cRVUqYQnZUTub3EjZQ1r449VR+MRnZKNuTsv4/9KV4QTEanFoSsp2Ho2ESKnybRn68tpU6rCQC2yjq5duxa7d++WyU7EHPXdxN5q0h0xfPTeM/UxfNExmdy+XwtfOX9NRKQGYjvW+79pkmENblULwZ7GkfJZ0aHvcePGYfDgwYiJiYG9vb3MV3r3jXSvU1338u1a/2N1LSJSEZHvQWzHqmZrgfFP1FG6OcbRoxZz0aLX3LNnT+23iCrtvafrY//lFOy7lIydF5LQlXNARKQwkea4rNb0293qsta0rnrUotdcuzaLQaiNn5sdXm7nL48/2HweeYVFSjeJiIzc7N/FItdC1PNyxIBQrmnSWaB+//33MX36dFniktRlbOdAuDtYIe5WDn44oNkGQUSkhLMJ6bLanzD92QYyQRPpKFD37dsXt2/flilDRXGOxx577J4bKcfeyhxTegbL4693XUFiei4vBxHpnFh0PP23cygpAZ5p7I1QfxdeBV3OUQ8dOhQnTpzAoEGD4OHhwWX2KvNckxpYdjgOJ6+mYdbWC5jTv6nSTSIiI/Nr+HUci70NGwszmdyEdByoN2/ejO3bt6Ndu3aP8NJUVcT+xOnPhuDZeQew4fR1DGxVCy38+GmWiHQjO6+wvFjQ6McD4O1sw1Ov66FvkTrU0ZH74NSsYU0n9G/hI4/fXX8WBUX37nUnIqoqX+68jMSMXPi62OKVDlx4rEig/uyzzzBp0iTExsY+cgOo6kx6MljuW4y8mYklh3itiKjqRSZmli9knd6rAawtmFZakUAt5qZFVrKAgAA4ODjAxcXlnhupQzU7S0wunRv6Yscl3EjnKn0iqtoFZO9uiEBRcQm6N/CUiZhIoTnqOXPmaOGlSRf6NPPBz8fi5cKyGZsuYN5Arsonoqrxy8mE8gVkU59hzQHFV31XxKxZszBy5Eg4OztX5mVIC0xNTTDjuYZ4eu5+bI64gX6XktGhTnWeWyLSqrScfMzcckEev9k1iAvItKhKi4GKutWpqalV+RJUAfW9HWWFLWHqxrPILWDGMiLSrk+3R+JWdj6C3O3xclvN+w3pQaAW8xWkDuOfCJIZy2Jv5eC7fdFKN4eIDMjp+DSsLM1A9sFzIbKiH2kPz6aRcLC2wLtPa+aMvt59BXG3spVuEhEZALFwTCwgE/2y55vWQKvarko3yeAwUBuRZxp5oW2gK/ILi/H+ryK1H0c8iOjRrDgSh7MJGXCwFumL6/F0VgEGaiPLWPa/XiGwMDPB7shkbIlIVLpJRKTHbmbk4tNtkfJ40pN1Ud3BSukmGSQGaiMTUN0eozoGyONpv55Dek6B0k0iIj01beM5ZOYVorGPM15qWUvp5hisKg3U7du3h40Nc7yqzehOgahd3Q4pWXmYtU2znYKI6GFsP5eIbecSYW5qglnPN2QJSzUG6qioKLz77rsYMGAAkpKS5GNbt27FuXPnyp+zZcsWeHl5aaelpDUipd/M3g3l8U9H4xEWfYtnl4gqLCO3QG71FF7tUBv1vFj7QXWBeu/evbIO9ZEjR7Bu3TpkZWXJx8PDwzFt2jRtt5GqQMvarhgQqina8d/1EdxbTUQVJualb2bkwc/VFm90CeKZU2Ognjx5MmbMmIEdO3bA0tKy/PHOnTsjLCxMm+2jKjS5Rz25+CM6ORvf7L7Cc01E/+pEXCqWH4mTxx/1bsiiG2oN1BEREejdu/ffHnd3d0dKSoo22kU64GRjgenPNpDH8/dG4dLNTJ53InqgvMIiTP5Fs2e6T7OaaBPoxrOl1kAtcnffuHHjb4+fOnUKNWrU0Ea7SEd6hHiiaz0PFBSVYPIvZ1BczL3VRHR/C/ZG43JSFlztLPF/T3HPtKoDdf/+/fHOO+8gMTFR7s0tLi7GwYMHMWHCBAwZMqTSjRFFPMT3GzduXKW/B1Vmb3UD2FmayQpbInkBEdFfXUnKwte7NFNkojKWs+2f056kwkAtim0EBwfDx8dHLiSrX78+OnTogDZt2siV4JVx7NgxLFiwAI0aNarUv6fK83a2waTumrrVH2+LREIa61YT0Z/ESNuUdWeQX1SMx+tWx7ONvXl61ByoRdpJ0ZP+6quvEB0djU2bNmH58uW4ePEili1bBjMzs4duhAj2AwcOxMKFC1GtWrWH/vf06Aa1qoXHfJ2RlVeIKevEHBSHwIlIY8nhWFlnWoy8zXguRI7EkcoDdWBgIK5duyZ71D179kTfvn0RFFT5JfpjxozBU089ha5du/7j8/Ly8pCRkXHPjbTDzNQEn/ZpLKve7LuUjDXHr/HUEhFiU7Lx8baL8kyIXN41q9nyrKg9UJuamsqgfOuWdpJkrFq1CidPnsTMmTP/9bniOU5OTuU38UGBtJtedEK3OvL4g03ncSOdQ+BExj7kPemXM8gtKEabAFe8FOqrdJOMkmllF31NnDgRZ89qMtNUVnx8PN58802sWLEC1tbW//r8KVOmID09vfwm/j1p14h2tdHU11nm7+UQOJFxWxYWh6MxqbC1NMPHLzSCqSmHvJVgUlKJyUgxj5yTk4PCwkKZ8OSv+bxTU1Mr9H02bNgg92PfPa9dVFQk5z9Ez10Mdf/TnLcY+hY9axG0HR2Zwk5briRloudXB2Q5zE9fbIQ+zTlyQWRsrt7KwZNz9uFOQRE+6NUAg1v7Kd0kg/Iw8cu8Mi8wZ84caEOXLl1k8pS7DR8+XK4oF9u/KrMwjR5doLsD3nqiDmZtvYj/bTqP9kHV4en07yMeRGRIQ97hMki3qu2CgayMpahKBeqhQ4dq5cUdHBwQEhJyz2N2dnZwdXX92+OkW6+0r41tZxNxOj5N5gL/YWhzrvQkMhIin0JYtGbI+9MXG3PIW9/LXObm5nIltgESq8Bn92kkV4HvupiEX04mKN0kItKB+NQczNyqWeU9uUcwfFy4ylsvA3V2djbGjh0rc3uLHrCYs7779ij27NmjtaF1evQh8PFdNavAp/92DteZCIXIoBUVl2DCmnDk5Behpb8LBnHIW38D9aRJk7Br1y7Mnz8fVlZW+P777zF9+nR4e3tj6dKl2m8lKeaV9v5o4uOMzNxC+QfMXOBEhuuHA9E4EpMqE5twyFvPA/Vvv/2Gb775Bi+88ALMzc3Rvn17mTpUpBYVW63IcJibmeKLfk1gY2GGQ1G38OPBGKWbRERV4MKNDMzefqk8l7evK4e89TpQi+1XtWvXlsdiWXnZdqx27dph37592m0hKc7fzQ7vPq2plPPJ9khEJrIcJpEhyS0owvifT8tc3qKaXl9uydT/QC2CdEyMpmcltlKtXr26vKctSmCS4REZiToHu8u91eN+Pi3r0hKRYfh8xyVcTMyEm70lZr3QkDs8DCFQi73O4eHh8njy5MmYN2+ezCw2fvx4mbGMDI9IQiP+gF3sLOUQ2Rc7LivdJCLSgsNRt7Bwf7Q8nvV8I7jZW/G8GkJmsr+Ki4vDiRMnZLEOXZapZGYy3RN7q0cuPwFRPOfnV1sj1N9FgVYQkTZk5Bag+xf7cD09FwNCfTDzeZYZVmP8euR91EKtWrXw/PPPs5a0Eege4ok+zWpCfLx7a/VpZOYWKN0kIqqk9zeek0G6lqst3n2qPs+jIWUm+9///veP/3/q1KmVbQ/pgWnPNkBYzC3Ep97B+7+ex2d9GyvdJCJ6SJvOXMe6UwkQdTY+79sEdlaVCgek1qHvpk2b3nO/oKBALi4TW7UCAgJk2Upd4NC3co7FpqLfgsMoLgG+7N8EvZrUULA1RPSw2cd6frVf5kd4vXMg3u5WlyfQ0IpynDp16r4vOmzYMFkNiwxfCz8XvN45CF/uvIz/W39WJkWp5WqndLOI6F8UFBXjjVWnZJAWJW3f6BLEc6ZyWpmjFsQnApGd7L333tPWtySVE5/EQ/1ckJVXiNd/OiW3bhGRun2x4xJOXU2Dg7U5vurfFBZmWgsDVEW0eoVEF17cyHiyls3p3wRONhY4cy0ds3+PVLpJRPQPDlxOwfy9UeVbsVhwQz9Uauj7q6++uue+mOa+ceMGli1bhh49emirbaQHvJ1t8OmLjfDqshP4bl80Wge4olNdd6WbRUR/kZKVh/GrT8sdGwNCffFUIy+eI0NeTObv73/PfVNTU1SvXh2dO3fGlClTZJ1pXeBiMvWYuvEslh6Og6udJba+2R7ujtZKN4mISoliOsMXH8PeS8mo42GPjWPawcbSjOfHkBeTlaUPJSrz3571cDQmVaYhfGt1OJa+HMpi80Qq8cOBGBmkrcxNMXfAYwzSeoarCEgrrC3M8PVLTWWVrQNX/pwHIyJlhcen4ZPtF+Xxe0/XR11P3Yx4kvZUqkcttmCJ3M8VsW7dusq8BOmhQHcHTH+2ASb9cgaf/R6Jx3yryTlrIlLG7ex8jF5xEgVFJejewBMDW/ryUhhLj1qMq+/cuRPHjx8vf0zk+t61a5ccaxf/v+xGxqVP85p44bGaMhGK2LKVlJGrdJOIjHZeWiweS0i7I1OEfvxiI1bFMqYetYeHB/r27Ytvv/0WZmaaBQlFRUUYPXq0DNSffvqptttJekKMtMx4LgTnrqfL+eoxK09i5SutuFeTSMfm7b6CPZGaeen5A5vJbZRkRD3qH3/8ERMmTCgP0oI4fuutt+T/I+MmVpPOH9QMDlbmOBZ7G59s08yPEZFu7L+cjM//uCSPP3guBPW9/3lVMRlgoC4sLMTFi39/8xWPFRczOxUB/m52+LSPpmTewv0x2Hb2Bk8LkQ5cT7uDN1dp9kv3a+6Dvs19eN6Nceh7+PDhGDFiBKKiohAaGiofO3LkCGbOnCn/H5HQPcQLr7T3l4F64pozqOvpKAM4EVUNkcZXTDelZuejvpcjpvdqwFNtrIF69uzZ8PT0xGeffSYzkgne3t6YNGkS3n77bW23kfTYpO7BOB2fJofARy0/gfWj23IPJ1EV+WjLhfI83t8Oaia3TZKRDn3n5eVh7NixSEhIQFpaGk6fPi0DdOPGje+ZtyYSCf+/fukxuNlbycVl/10fIVPOEpF2bTydgMWHYuWxqC/t62rLU2zMgbpXr15YunSpPBZz0t26dcPnn3+O5557DvPnz9d2G0nPeThaY+6ApjAzNcH6UwkySxIRac/ZhHRMWntGHo96PABP1Pfg6TX2QH3y5Em0b99eHq9du1Zu14qLi5PB+68FO4gEkfjk3afqlQ/PiVWpRKSdYhuvLj2OvMJiPF63OiZ0q8vTamAqFahzcnLKC2/8/vvveP7552VhjlatWsmATXQ/w9r4oU8zTTKUsStPIe5WNk8U0SMuHhu9/CSup+eitpsdvuyvGbkiw1KpQB0YGIgNGzYgPj4e27dvl0PfQlJS0r9WASEjT4bSOwRNfJyRfqcAryw9jqy8QqWbRaS3/rfpHI7GpsLeyhzfDWFSE0NVqUA9depUmfDEz88PLVu2ROvWrct7102bNtV2G8mAWJmbYcHgZnB3sMKlm1l46+fTMtUhET2cn45exfKwqxBlF+b0ayJz7ZNhqlQ9aiExMVFuzRIrvcWwt3D06FHZow4ODoYusB61/jp59Tb6LwhDflEx3uwShPFP1FG6SUR643hsKgYsDJPFNiZ0q4OxnYOUbhJVYfyqdJlLsY9a9J7LgrQgkp/oKkiTfhOVtUROcOHLnZfl1hIi+ndibcery07IIN2zoSfGdArkaTNwrEdNiunbwkdmLhNE5jLRSyCiB0vPKcDLi4/JzGMNazhhdp/GrIhlBBioSVGTe9STez7FELjoJVy9lcMrQnQfBUXFGLXiBKKSs+HlZI3vhzaHrWWlkkuSnmGgJkWJrSRf9m+CkBqOspcwfPFRuSKciP4klhK9u/4sDkXdgp2lGX4Y2kImEiLjwEBNihO9AvHG4+loLXsLo1eI+TdWYSMqs2BfNH4+Hg+xRXruS01ZttLIMFCTKojewQ/DxFCeGQ5euYX3NpxlTnAiAFsjbmDWVk1Z4alP10fnYKYHNTYM1KQaDbydZE5w0WtYdSwec3ddUbpJRIo6FpuKN38+LY+Htq6FYW01iy/JuDBQk6p0qeeB6c9qauh+vuOSTOpAZIwu3czEiMXHZJrQrvU88N7T9ZVuEimEgZpUZ3BrP4wt3Rv6f+sjsOP8TaWbRKRT19PuYOiPR5GRW4hmtarJkSZzM75dGytFr/zMmTPRokULWeDD3d1dlsmMjIxUskmkEm93q4O+zcsKeJzEiTjusSbjkJaTL4P0jfRcBLrb44ehzWFjaaZ0s8hYA/XevXsxZswYhIWFYceOHSgoKJAFPrKzWVXJ2IkCHh/1bojOwe6yfN/Li4/j8s1MpZtFVKVyC4rwnyXHcTkpCx6OVljyciicbS151o1cpXN9V4Xk5GTZsxYBvEOHDv/6fOb6Nnw5+YV4aeERnI5Pg7eTNdaOagNvZxulm0WkdYUyoclJOdXjYG2ONSNbI9iT1QgNlU5yfVcF0WDBxcXlvv8/Ly9P/nB338jw91j/OKwFale3kzV3B31/BMmZeUo3i0irRAW5iWvPyCBtaW6KhUOaM0iT+gJ1cXExxo0bh7Zt2yIkRFOs4X5z2uITSNnNx8dH5+0k3XOxs8SyES1Rw9kG0SnZGPzDETmPR2QwWcc2nsX6UwkwNzXBNy89hla1XZVuFqmIaoa+R40aha1bt+LAgQOoWbPmA3vU4lZG9KhFsK7I0AHpv9iUbPRZcFj2qBv7OGPFf1rC3oq5jkl/ibffj7ZcwML9MbKu9Jf9m+LZxt5KN4t0QO+GvseOHYtNmzZh9+7dDwzSgpWVlfyB7r6R8fBzs8PyES1RzdYC4fFpsorQnfwipZtFVGmixKsI0sLHzzdikCb1BWrxaVIE6fXr12PXrl3w92fWHfpndT0dsPTllnCwMsfRmFSMXH4CeYUM1qR/Fu6Lxpw/Lsvjac/Ul2VfiVQXqMXWrOXLl2PlypVyL3ViYqK83blzR8lmkco1rOmEH4e3gI2FGfZeSsaYFScZrEmvLDoYgw+3XJDHE5+si+FMDUpqnaMWe2XvZ9GiRRg2bNi//ntuzzJuBy6nYMSSY3KfdZdgd3wz6DFYmTMxBKnb9/ujMWOzJkiP6RSAiU8GK90kUoDezFGLzwj3u1UkSBO1C3KT5TGtzE2x82ISRi47IRNGEOlDkBZpcid0q6t0k0gPqGIxGdGjBGuxz9rawhS7I5PxGoM1qdSCvVHlQfqNzoEyTe6DRhWJ7sZATXqvbaAmWJfNWb+y9Dh71qQq8/dEYWZpTelxXYPwVre6DNJUYQzUZBDaBLhhUekCs/2XU+TWray8QqWbRUZOTOV9+cdlfLxNE6THd62DcV3rKN0s0jMM1GQwRDanxcNbwM7SDIeibmHgwjDczmYGM1IuLej0387jiz8uyfsTutXBm12DeDnooTFQk0FpWdsVK19ppUmKci0dfRccRmJ6rtLNIiNTUFSMCWvCsfhQrLz//jP1MbYzgzRVDgM1GRyRXnT1a63h6WgtywW+MP+QTD9KpAti58Go5Sex7lQCzExN8EW/xhjGfdL0CBioySAFeTjIMoF+rrZISLuDF789jPPXWW2NqlZmbgGGLTqKPy5oqmAtGNQMvZs+OC0yUUUwUJPB8nGxxZqRbVDPyxEpWXnot+AwDl5JUbpZZKDEFEu/BWEIi06VxWKWvhyKrvU9lG4WGQAGajJo1R2ssOrVVgj1c0FmXiGG/ngUa09cU7pZZGAuJmag9zcHcf5GBtzsLfHTK61YqpK0hoGaDJ6TjQWWjgjFM429UVhcIhf5fLHjktw6Q/So9l9OxovzD+NGei5qV7fD+tFtZT56Im1hoCajYG1hhi/7NcHoxwPKywu+vSYc+YXFSjeN9NjqY/EYvkizZ7+lvwvWjWojp1yItImBmoyGqakJJnUPxsznG8rVuOtOJsih8LQc7rWmh98jPXt7JCb9ckaO0jzXxFuO2jjbWvJUktYxUJPRGRDqK1OOisQoh6Nv4dmvD8o5RqKKrux+ddlxfL37irz/eudAfNGvCSu3UZVhoCaj1LFOdfwyWgxT2uBqag6e/+YQtp29oXSzSOWik7PQ+5tD+ONCktx+9VmfxnibebupijFQk9EK9nTEr2PaoW2gK3LyizBy+Ul8/nukHNYk+qvdkUnoNe8griRlyWQ6a15rjReacY80VT0GajJq1ewssWR4KEa085f3v9p1RQ5rZuQWKN00UgmxO+CbPVdkoZfM3EI0r1UNv77eVmbAI9IFBmoyeuZmpnjv6fpyGFMMZ4phzae/OoAz19KM/twYO7HQUJRN/WRbJMRuPrG+QeSSd3ewVrppZEQYqIlKiWHMtSNbo2Y1zby1yBG+6GAM91sbqRNxt/HUVwfK56M/6t1Q7hgQx0S6xN84ors0qumMzW+0x5MNPFBQpClTOHL5CaTncCjcmIa6F+6LlilnRZ54kS9+/eg2eKmlr9JNIyNlUqLH6ZkyMjLg5OSE9PR0ODo6Kt0cMiDiz2LJoVh8tOUi8ouKZS/7qwFN8ZhvNaWbRlUoKTMX76w9g92RyfL+0428ZC/awdqC550Ui18M1ET/IOJaOsasPCmHwk1NgDGdAvF65yAOfxqg7ecSMWVdBFKz8+X1nfp0fQxs6QsTExOlm0YGiIGaSJt/ULkFmLrhLDacvi7vh9RwxBd9m8hSmqT/RPrP6b+ew5rSYi2i2tqcfk1Q15PXl6oOAzVRFdh85gb+b0ME0nIKZI9r0pN1Mbytv0xHSvrp0JUUvLPuDOJT70B0nF/rEIDxTwQxyxhVOQZqoiqSlJGLd375cw6ziY8zZr3QUCZPIf0hFgd+uOU8Vh/X9KJrONvINKCh/i5KN42MRAbnqImqdqHZqmPx+GjzBVnj2tzUBCM7BmBs50BZpYvUbWvEDUz99RySM/Pk/SGta2Hik3W5YIx0ioGaSAcS03Mx7dez2H7uprzv72aHD58LQZtAN55/FYpPzZHb7f64oLleAdXt8PELjdDcj71o0j0GaiIdEsU8pm48h6TSHlrPhp74b896qFmNdYnV4E5+EebvjcK3e6Nk/XExAjLq8QC5gp8jIKQUBmoiXf/R5RbI+sTLw+IganpYmZvitY4BGNUxADaWHA5Xaopi29lEzNh8QSYuEdoEuOL9ZxugDlfsk8IYqIkUcuFGBqb/dg5h0anyvreTtSyD+FzTGlwdrkMn4lIxa+tFHIu9XX4d3n26PnqEeHJfNKkCAzWRwj25rWcT8eFdPbk6HvaY0K0unqjvwUBRhS7fzMQn2yOx47xmHlqMbLzaoTZGPx7IkQ1SFQZqIhXILSiSaUi/2ROF9DuaXOFNfZ0xsVtdtA5wZcDWopiUbMzbfQXrTl6TUw9ia3vf5j54s2sQvJxstPlSRFrBQE2kIiJIf7cvCj8eiMWdgqLygD3m8UB0qefOgP0ILt3MlAH6t/DrMkALoqCK2G4V6M7MYqReDNREKi34MG/XFfx0LF6uPhaCPR3kCuSnGnrJuthUMaeu3saCvdHYdi6x/LHOwe5yLzsLp5A+YKAmUnnAFr1rsUJc5JkWvJysMahVLfRr4QM3eyulm6hK4sPNlogbWHQoFuHxaeWPiwViYqtVSA0nRdtH9DAYqIn0ZEh82eFYLDoYi1vZ+fIxSzNTPN3YC4Nb1ZLpSVm5Cbh2Owdrjl/DyqNXy7OJifP0TGNvvNaxNrdakV5ioCbSI3mFRbLgx5LDcff0FAPd7fH8YzXQu2kNo1sQJZKUbDt3QwboQ1G3yh/3cLSSH2L6h/py5IH0GgM1kZ46HZ+GpYdiseXsDeQWaOaxRVWntgFueKaxF7rW84CrgQ6Ni1Xy+y+nYOvZG/j93M3yaYGyRCUiOIthbgvO5ZMBYKAm0nOZuQXYGpGItSev4WiMJnmKILYdidzUTzbwRLf6HvBxsdX7KlYHrmiC8+6LScjO16yKF3xcbPDiYz54oVkNpmMlg8NATWRgxSQ2nk6QK5zPJmTc8//8XG1lEZB2gW5oXdsV1ewsofYFYWLFtgjOovd85lpa+baqskV14kNIz4ZeaF6rGkxZ65sMlN4F6nnz5uHTTz9FYmIiGjdujLlz5yI0NFSrPyiRoSysEsPC288l4lhs6j1BTgyR13F3QGMfJzT2cUbjms6o6+mg2FCxeGsRmdnEcP6pq+J2G2evZ5RvTbt7Lr5LsDu6h3jKNjM4kzHI0KdA/fPPP2PIkCH49ttv0bJlS8yZMwdr1qxBZGQk3N3d//HfMlCTsRcCORKdioNXUnAoKgWXbmb97TmW5qao7WaHgOr2sqxjgLs9fF1s4elkLRdjPWoQF28fqdn5uJGeK28xKVm4fDMLl5OycCUp65555jIudpZyBKBdkBvaB7kZ3UI5Ir0L1CI4t2jRAl9//bW8X1xcDB8fH7z++uuYPHnyP/5bBmqie/dnh8eny5Xj4dfS5NeM3L8Hyrt74K52VnIltaO1BeyszGBnZS5v1uZmKBH/lWiCsXiTyMkvknPnIvhm5hYiLacAiRm5f+sh302UlKzn5Si3molsbE19q8nhem47I2OX8RCB2hwKys/Px4kTJzBlypTyx0xNTdG1a1ccPnz4b8/Py8uTt7t/UCLScHewxhP1xc1D3i8uLsG123cQlZxVfhO93ITbd2Tt7MLiEqRk5cnboxK9czG/LHrrQR72CHJ3kF/9XO1kr56IKk/RQJ2SkoKioiJ4eGjeWMqI+xcvXvzb82fOnInp06frsIVE+kvM9fq62spbp+B7p5FEEE/NyUdieq5MIpKZV4js8lsRcguLYCJXmZvInrc4trY0g4O1BRyszOFgbQ5HGwt4OlrDw9GawZjIUAP1wxI977feeuueHrUYJieihw/iohfMdKVE6qdooHZzc4OZmRlu3tTUji0j7nt6ev7t+VZWVvJGRERkLBSdPLK0tESzZs2wc+fO8sfEYjJxv3Xr1ko2jYiISBUUH/oWQ9lDhw5F8+bN5d5psT0rOzsbw4cPV7ppREREilM8UPfr1w/JycmYOnWqTHjSpEkTbNu27W8LzIiIiIyR4vuoHwX3URMRkaHHL25wJCIiUjEGaiIiIhVjoCYiIlIxxReTPYqy6XWmEiUiIn1SFrcqskxMrwN1Zmam/MrsZEREpK9xTCwqM9hV3yI5yvXr1+Hg4KCVajxlKUnj4+ONtr41zwHPA38X+DfB94Wqf28UoVcEaW9vb1mMymB71OKHq1mzpta/r7gIxhqoy/Ac8Dzwd4F/E3xfqNr3xn/rSZfhYjIiIiIVY6AmIiJSMQbqu4jKXNOmTTPqCl08BzwP/F3g3wTfF9T13qjXi8mIiIgMHXvUREREKsZATUREpGIM1ERERCrGQE1ERKRiDNQP8Oyzz8LX1xfW1tbw8vLC4MGDZRY0YxIbG4sRI0bA398fNjY2CAgIkKse8/PzYUw+/PBDtGnTBra2tnB2doaxmDdvHvz8/OTfQMuWLXH06FEYk3379uGZZ56RmaNE5sMNGzbA2MycORMtWrSQ2R/d3d3x3HPPITIyEsZk/vz5aNSoUXmik9atW2Pr1q06bQMD9QN06tQJq1evlr+Uv/zyC6KiovDiiy/CmFy8eFGmaV2wYAHOnTuHL774At9++y3++9//wpiIDyZ9+vTBqFGjYCx+/vlnvPXWW/KD2cmTJ9G4cWM8+eSTSEpKgrHIzs6WP7f4wGKs9u7dizFjxiAsLAw7duxAQUEBunXrJs+NsahZsyZmzZqFEydO4Pjx4+jcuTN69eol3xN1RmzPon+3cePGEhMTk5L8/HyjPl2ffPJJib+/f4kxWrRoUYmTk1OJMQgNDS0ZM2ZM+f2ioqISb2/vkpkzZ5YYI/FWuX79+hJjl5SUJM/F3r17S4xZtWrVSr7//nudvR571BWQmpqKFStWyOFPCwsLGLP09HS4uLgo3Qyq4hEE0Xvo2rXrPXn1xf3Dhw/z3Bv5379grO8BRUVFWLVqlRxREEPgusJA/Q/eeecd2NnZwdXVFVevXsXGjRthzK5cuYK5c+fitddeU7opVIVSUlLkG5KHh8c9j4v7iYmJPPdGSkyDjRs3Dm3btkVISAiMSUREBOzt7WVWspEjR2L9+vWoX7++zl7fqAL15MmT5aKQf7qJedkyEydOxKlTp/D777/DzMwMQ4YMqVCRb0M7D0JCQgK6d+8u52pfeeUVGOM5IDJmYq767NmzskdpbOrWrYvTp0/jyJEjcq3K0KFDcf78eZ29vlGlEE1OTsatW7f+8Tm1a9eGpaXl3x6/du2arEV66NAhnQ55qOE8iNXujz/+OFq1aoXFixf/a+1UQ/1dED+76FGkpaXB0Ie+xQr3tWvXylW+ZcSbk/jZjXFkSXxwE72ou8+HMRk7dqy87mIlvNgFYuy6du0qd8GIhba6oNf1qB9W9erV5a2ywz5CXl4ejOk8iJ60WAHfrFkzLFq0yCCC9KP+Lhg68eFEXO+dO3eWBybx+y/uizdsMh6iH/f666/LDyl79uxhkC4l/h50GQuMKlBXlBjeOHbsGNq1a4dq1arJrVnvvfee/ASl773phyGCtOhJ16pVC7Nnz5a90DKenp4wFmJ9glhQKL6KuVsxBCYEBgbKeStDJLZmiR508+bNERoaijlz5sgFNMOHD4exyMrKkusyysTExMhrLxZSiRwLxjLcvXLlStmbFnupy9YoODk5ydwKxmDKlCno0aOHvOaZmZnyfIgPLdu3b9ddI3S2vlyPnDlzpqRTp04lLi4uJVZWViV+fn4lI0eOLLl27VqJsW1HEr8i97sZk6FDh973HOzevbvEkM2dO7fE19e3xNLSUm7XCgsLKzEm4vre77qL3wdj8aC/f/HeYCxefvnlklq1asm/g+rVq5d06dKl5Pfff9dpG4xqjpqIiEjfGMaEIxERkYFioCYiIlIxBmoiIiIVY6AmIiJSMQZqIiIiFWOgJiIiUjEGaiIiIhVjoCYiIlIxBmoiIiIVY6AmIiJSMQZqIiIiFWOgJqJyokKaqIz20UcflT8marCL0peizCUR6R6LchDRPbZs2SLrUIsAXbduXTRp0gS9evXC559/zjNFpAAGaiK6bx3iP/74Q9ajjoiIkPXZrayseKaIFMBATUR/c+fOHYSEhCA+Ph4nTpxAw4YNeZaIFMI5aiL6m6ioKFy/fh3FxcWIjY3lGSJSEHvURHSP/Px8hIaGyrlpMUc9Z84cOfzt7u7OM0WkAAZqIrrHxIkTsXbtWoSHh8Pe3h4dO3aEk5MTNm3axDNFpAAOfRNRuT179sge9LJly+Do6AhTU1N5vH//fsyfP59nikgB7FETERGpGHvUREREKsZATUREpGIM1ERERCrGQE1ERKRiDNREREQqxkBNRESkYgzUREREKsZATUREpGIM1ERERCrGQE1ERKRiDNRERERQr/8HH+hTL9RLV9QAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def square_me(x):\n", " return x**2\n", "\n", "\n", "x = np.arange(-3, 3.01, 0.01)\n", "\n", "plt.figure(figsize=(5, 3.2))\n", "plt.plot(x, square_me(x))\n", "plt.xlabel(\"x\"); plt.ylabel(\"square_me(x)\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3509139f", "metadata": {}, "source": [ "### Writing our own average\n", "\n", "R returns the last expression automatically; **Python needs an explicit `return`**." ] }, { "cell_type": "code", "execution_count": 35, "id": "1cc5bb07", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "50.5 50.5\n", "True\n" ] } ], "source": [ "def avg(x):\n", " s = sum(x)\n", " n = len(x)\n", " return s / n # R would return s/n implicitly -- Python will not\n", "\n", "\n", "x = list(range(1, 101))\n", "print(avg(x), np.mean(x))\n", "print(np.isclose(np.mean(x), avg(x))) # R: identical(mean(x), avg(x))" ] }, { "cell_type": "markdown", "id": "7775eaf9", "metadata": {}, "source": [ "### Variables inside a function are local" ] }, { "cell_type": "code", "execution_count": 36, "id": "48ddc966", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.5\n", "3\n" ] } ], "source": [ "s = 3\n", "print(avg(range(1, 11))) # the s inside avg() is a different s\n", "print(s) # ... the global s is untouched" ] }, { "cell_type": "markdown", "id": "80bed105", "metadata": {}, "source": [ "### Default arguments\n", "\n", "R:\n", "```r\n", "avg <- function(x, arithmetic = TRUE){\n", " n <- length(x)\n", " ifelse(arithmetic, sum(x)/n, prod(x)^(1/n))\n", "}\n", "```" ] }, { "cell_type": "code", "execution_count": 37, "id": "c80a7c6d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "arithmetic mean: 2.5\n", "geometric mean: 2.213363839400643\n" ] } ], "source": [ "def avg(x, arithmetic=True):\n", " x = np.asarray(x, dtype=float)\n", " n = len(x)\n", " return x.sum() / n if arithmetic else np.prod(x) ** (1 / n)\n", "\n", "\n", "print(\"arithmetic mean:\", avg([1, 2, 3, 4]))\n", "print(\"geometric mean: \", avg([1, 2, 3, 4], arithmetic=False))" ] }, { "cell_type": "code", "execution_count": 38, "id": "158fc950", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5050" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def compute_s_n(n):\n", " return sum(range(1, n + 1))\n", "\n", "\n", "compute_s_n(100)" ] }, { "cell_type": "code", "execution_count": 39, "id": "89eb1588", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE1CAYAAAAoBaPAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAI9BJREFUeJzt3XtsVGX6wPGnLW2hXGTLpZelIFrlDhoQaLisQKVcJCKYiLJatCmBpUSoIKmRmxrLgvG6KH+sgpsIKhvRiBHtltuCBd0iUVC7QoiUpRfUFLlIW9r55Xl15teR6b0zc86c7yc5mZ45p9Pp4dBn3vd93ucNc7lcLgEAAJYUHuw3AAAA6kegBgDAwgjUAABYGIEaAAALI1ADAGBhBGoAACyMQA0AgIW1C/YbsILa2lo5e/asdO7cWcLCwoL9dgAAIc7lcsmFCxckMTFRwsMbbjMTqEVMkE5KSgrUvw8AAEZxcbH06tVLGkKgFjEtafcF69KlS4MXDACA1vr5559NA9EdfxpCoBbxdHdrkCZQAwACpSnDrSSTAQBgYQRqAAAsjEANAICFEagBALAwAjUAAE1wpbpGin+6bB4DiaxvAAAaUFPrkufz/iuvHTglv1TXSIfICMkY21eW3nGzRIT7v0gWgRoAgAZokP7bnhOefQ3W7v1laf3E3+j6BgCgHtrNrS1pX14/eCog3eAEagAA6nHuQqVpQftyuarGHPc3AjUAAPXo0TnajEn7EhMVYY77G4EaAIB6tP8tccyXh8f0Ncf9jWQyAAAaoNnd7jFp7e7WlrQGaffz/hbm0kUxHU5XMbnuuuvk/PnzLMoBAPBJE8d0TFq7u1vbkm5O3KFFDQBAE2hwToqNkUBjjBoAAAsjUAMAYGEEagAALIxADQCAhRGoAQCwMAI1AAAWRqAGAMDCCNQAAFhYUAP1q6++KkOHDjVVWXRLSUmRjz76yHP8ypUrsmjRIunWrZt06tRJZs+eLWVlZV6vcfr0aZk+fbrExMRIz549Zfny5XL16tUg/DYAAKtXFiv+6XJAlqZsS0GtTNarVy9Zt26d3HTTTaKVTN944w2566675IsvvpBBgwbJ0qVL5cMPP5Tt27ebUmtZWVkya9YsOXjwoPn+mpoaE6Tj4+Pl008/lZKSEnnwwQclMjJSnnnmmWD+agAAi6ipdcnzef8160rrkpUdfltoQ2t1R4SHidVZrtZ3bGysbNiwQe655x7p0aOHbN261Xytvv32WxkwYIAUFBTI6NGjTev7zjvvlLNnz0pcXJw5Z9OmTbJixQo5d+6cREVFNelnUusbAELXsx8Xyd/2nLjm+awJybIsrV9Q3lNz4o5lxqi1dfzWW2/JpUuXTBd4YWGhVFdXS2pqquec/v37S+/evU2gVvo4ZMgQT5BWaWlp5gIcP3683p9VWVlpzqm7AQBCz5XqGtOS9kVXw7JDN3jQA/VXX31lxp+jo6NlwYIFsmPHDhk4cKCUlpaaFnHXrl29ztegrMeUPtYN0u7j7mP1yc3NNZ9k3FtSUpJffjcAQHCdu1Bpurt90SUr9bjVBT1Q9+vXT44ePSqHDx+WhQsXSnp6unz99dd+/Zk5OTmmu8G9FRcX+/XnAQCCo0fnaDMm7YuuK63HrS7ogVpbzcnJyTJ8+HDT0h02bJi8+OKLJkGsqqpKKioqvM7XrG89pvTx91ng7n33Ob5o692dae7eAAChuTRlxti+Po89PKZvq9eVdkSg/r3a2lozhqyBW7O38/PzPceKiorMdCwdw1b6qF3n5eXlnnPy8vJM4NXucwAAlt5xs0kc0xa00kfd1+ftIKjTs7QLeurUqSZB7MKFCybDe+/evfLxxx+bseOMjAzJzs42meAafBcvXmyCs2Z8q8mTJ5uA/MADD8j69evNuPQTTzxh5l5rqxkAgIjwMJPdnTUx2YxJa3e3HVrSlgjU2hLWec86/1kDsxY/0SB9xx13mOPPP/+8hIeHm0In2srWjO5XXnnF8/0RERGyc+dOM7atAbxjx45mjPvJJ58M4m8FALCi9pERkhQbI3ZjuXnUwcA8agBAINlyHjUAALgWgRoAAAsjUAMAYGEEagAALIxADQCAhRGoAQC2ccWma0rbdh41AABOWFO6NQjUAADLez7vv15rSmuwdu8Ha03pQKHrGwBgaVdCYE3p1iBQAwAs7VwIrCndGgRqAICl9QiBNaVbg0ANALC09iGwpnRrkEwGALC8pb+tHa1j0trdrS1pDdJ2WVO6NVg9i9WzAMA2rlTX2HJN6dasnkWLGgBgG+1tuqZ0azBGDQCAhRGoAQCwMAI1AAAWRqAGAMDCCNQAAFgYgRoAAAsjUAMAAsaJ60m3FvOoAQB+5+T1pFuLQA0A8Dsnrydt667v3Nxcue2226Rz587Ss2dPmTlzphQVFXmdc/vtt0tYWJjXtmDBAq9zTp8+LdOnT5eYmBjzOsuXL5erV68G+LcBAPji9PWkbR2o9+3bJ4sWLZJDhw5JXl6eVFdXy+TJk+XSpUte52VmZkpJSYlnW79+vedYTU2NCdJVVVXy6aefyhtvvCFbtmyRVatWBeE3AgD8ntPXk7Z11/euXbu89jXAaou4sLBQxo8f73leW8rx8fE+X+OTTz6Rr7/+Wv71r39JXFyc3HLLLfLUU0/JihUrZM2aNRIVFXXN91RWVpqtbnF0AIB/15P2FaydsJ50SGV96yoiKjY21uv5N998U7p37y6DBw+WnJwcuXz5sudYQUGBDBkyxARpt7S0NBN8jx8/Xm+Xu65a4t6SkpL89jsBgNM5fT3pkEkmq62tlSVLlsiYMWNMQHa7//77pU+fPpKYmChffvmlaSnrOPa7775rjpeWlnoFaeXe12O+aLDPzs727GtQJ1gDgP84eT3pkAnUOlZ97NgxOXDggNfz8+fP93ytLeeEhASZNGmSnDx5Um688cYW/azo6GizAQACQ6dgaXZ31sTkkFhP2nFd31lZWbJz507Zs2eP9OrVq8FzR40aZR5PnPg1rV/HrsvKyrzOce/XN64NAAjuetIEaZsEapfLZYL0jh07ZPfu3dK3r+8xjLqOHj1qHrVlrVJSUuSrr76S8vJyzzmaQd6lSxcZOHCgH989AAAh3vWt3d1bt26V999/38yldo8pa4JXhw4dTPe2Hp82bZp069bNjFEvXbrUZIQPHTrUnKvTuTQgP/DAA2balr7GE088YV6b7m0AgN2FubRZG6wfHua7bNzmzZtl3rx5UlxcLH/+85/N2LXOrdaEr7vvvtsEYm0xu33//feycOFC2bt3r3Ts2FHS09Nl3bp10q5d0z6HaDKZfjjQrPO6rwsAgD80J+4ENVBbBYEaAGDVuGOJZDIAAOAbgRoA0GQsU+ngedQAAOtimcrgIVADABrFMpXBQ9c3AKBBLFMZXARqAECDWKYyuAjUAIAmLVPpC8tU+h+BGgDQIJapDC6SyQAAjWKZyuChMhmVyQCgWYllLFMZ2MpktKgBAM1ephKBwxg1AAAWRqAGAMDCCNQAAFgYgRoAAAsjUAMAYGEEagBwGJaqtBemZwGAQ7BUpT0RqAHAIViq0p7o+gYAB2CpSvsiUAOAA7BUpX0RqAHAAViq0r4I1ADgACxVaV9BDdS5ubly2223SefOnaVnz54yc+ZMKSoq8jrnypUrsmjRIunWrZt06tRJZs+eLWVlZV7nnD59WqZPny4xMTHmdZYvXy5Xr14N8G8DANZfqjJrQrLEREWYfX3UffcSlrCmoC5zOWXKFJkzZ44J1hpYH3/8cTl27Jh8/fXX0rFjR3POwoUL5cMPP5QtW7aYJcGysrIkPDxcDh48aI7X1NTILbfcIvHx8bJhwwYpKSmRBx98UDIzM+WZZ55p8+XGAMDuWKoy+JoTdyy1HvW5c+dMi3jfvn0yfvx48wv06NFDtm7dKvfcc48559tvv5UBAwZIQUGBjB49Wj766CO588475ezZsxIXF2fO2bRpk6xYscK8XlRUVKM/l0ANAAik5sQdS41R6xtWsbGx5rGwsFCqq6slNTXVc07//v2ld+/eJlArfRwyZIgnSKu0tDRzEY4fP+7z51RWVprjdTcAAKzIMoG6trZWlixZImPGjJHBgweb50pLS02LuGvXrl7nalDWY+5z6gZp93H3sfrGxvWTjHtLSkry028FAECIBGpNGNPx6bfeesvvPysnJ8e03t1bcXGx338mAAC2LSGqCWI7d+6U/fv3S69evTzPa4JYVVWVVFRUeLWqNetbj7nP+eyzz7xez50V7j7n96Kjo80GAIDVBbVFrXlsGqR37Nghu3fvlr59+3odHz58uERGRkp+fr7nOZ2+pdOxUlJSzL4+fvXVV1JeXu45Jy8vzwzODxw4MIC/DQAEBqtfOUu7YHd3a0b3+++/b+ZSu8eUddy4Q4cO5jEjI0Oys7NNgpkG38WLF5vgrBnfavLkySYgP/DAA7J+/XrzGk888YR5bVrNAEIJq185U1CnZ4WFhfl8fvPmzTJv3jxPwZNHH31Utm3bZrK1NaP7lVde8erW/v777818671795r51+np6bJu3Tpp165pn0OYngXADp79uEj+tufENc9r0ZJlaf2C8p7gsHnUwUKgBmCH7u5bn8yTX6prrjmmFcaOrLzDlAmFPdh2HjUAwDdWv3IuAjUA2ACrXzkXgRoAbIDVr5zLEvOoAQCNc69y9frBU3K5qsaMTT88pi+rX4U4kslIJgNgM6x+5axkMlrUAGDDbvCk2Jhgvw0ECGPUAABYGIEaAAALI1ADAGBhBGoACDAW1UBzkEwGAAHCohoIWKCuqamRLVu2mOUndXnJ2tpar+O6ZCUAwNvzef/1WlRD63a791lUA20aqB955BETqKdPny6DBw+udxUsAMD/d3e/duCUz8uhBUyyJiazqAbaLlC/9dZb8s4778i0adNa8u0A4DhNWVSDudFos2SyqKgoSU5Obsm3AoAjsagGAhqoH330UXnxxReFpawBoGlYVAMB7fo+cOCA7NmzRz766CMZNGiQREZGeh1/9913W/yGACBUsagGAhaou3btKnfffXeLfiAAOFVEeJjJ7tbEMR2T1u5wbWkDQVs96+DBgzJixAiJjo6WUFnFBACAQMYdv1Ymmzp1qvzvf//z548AACCk+TVQk2wGAEDrUOsbAFqAet0IFGp9A0AzUK8bgUagBoBmoF43Qqrru7Ea4Pv375cZM2ZIYmKiOfe9997zOj5v3jzzfN1typQpXuf89NNPMnfuXJM1p9PGMjIy5OLFi375fQA4W2P1uvU4YIlA/csvv8jly5c9+99//7288MIL8sknnzQrmezSpUsybNgw2bhxY73naGAuKSnxbNu2bfM6rkH6+PHjkpeXJzt37jTBf/78+S35tQCg1fW6AUt0fd91110ya9YsWbBggVRUVMioUaNMdbIffvhBnnvuOVm4cKE578KFC41O39KtIToHOz4+3uexb775Rnbt2iWff/65ma+tXn75ZbNYyLPPPmta6r5UVlaare58NgBoar1uX8E6JirCHAcs0aI+cuSIjBs3znz9z3/+U+Li4kyr+h//+Ie89NJLbfoG9+7dKz179pR+/fqZDwA//vij51hBQYHp7nYHaZWamirh4eFy+PDhel8zNzfXTDR3b0lJSW36ngGEJup1wzaBWru9O3fubL7W7m5tXWtwHD16tAnYbUW7vTX45+fny1//+lfZt2+faYHX1Pz6aba0tNQE8bratWsnsbGx5lh9cnJyTDUY91ZcXNxm7xlA6NfrzpqQbFrQSh91313HG7BE17cucamJX1rv++OPP5alS5ea58vLy9u0BOecOXM8Xw8ZMkSGDh0qN954o2llT5o0qcWvq93pVi9rCsCaqNcNW7SoV61aJcuWLZPrr7/ejE+npKR4Wte33nqr+MsNN9wg3bt3lxMnTph9HbvWDwd1Xb161WSC1zeuDQBt1Q2eFBvDohqwZov6nnvukbFjx5osbM3adtNWrj9X1Tpz5owZo05ISDD7+gFBk9kKCwtl+PDh5rndu3dLbW2t+QABAIBjC55oi/X3rdaRI0c26zV0vrO7daxOnTolR48eNWPMuq1du1Zmz55tfs7JkyflscceM93uaWlp5vwBAwaYcezMzEzZtGmTVFdXS1ZWlukyry/jGwAAO/HrMpeN0bHmCRMmXPN8enq6vPrqqzJz5kz54osvTKtZA+/kyZPlqaeeMlnmbtrNrcH5gw8+MAltGtg187xTp05Nfh8scwk4jxYnYU1oBEtz4k5QA7VVEKgB56BWN+wWd6j1DcBRqNUNu2GZSwCOQa1u2BGBGoBjUKsbdkSgBuC4Wt2+UKsbVkWgBuAY1OqGHZFMBsBR3DW5df1oXZpSW9IPj+lLrW5YFtOzmJ4FOBLzqBFMTM8CgCbW6gasjjFqAAAsjEANwLZd18U/XTaPQCgjmQyArVACFE5DoAZgK5QAhdPQ9Q3ANigBCiciUAOwDUqAwokI1ABsgxKgcCICNQDboAQonIhkMgC2QglQOA0lRCkhCtgSJUBhZ5QQBRDyKAEKp2CMGgAACyNQAwgayoACjSOZDEDAUQYUaDoCNYCAowwoYJOu7/3798uMGTMkMTFRwsLC5L333vM67nK5ZNWqVZKQkCAdOnSQ1NRU+e6777zO+emnn2Tu3LnSpUsX6dq1q2RkZMjFixcD/JsAaCrKgAI2CtSXLl2SYcOGycaNG30eX79+vbz00kuyadMmOXz4sHTs2FHS0tLkypUrnnM0SB8/flzy8vJk586dJvjPnz8/gL8FgOagDChg03nU2qLesWOHzJw50+zr29KW9qOPPirLli0zz50/f17i4uJky5YtMmfOHPnmm29k4MCB8vnnn8uIESPMObt27ZJp06bJmTNnzPe39Xw2AK1vUd/6ZJ784mMd6ZioCDmy8g4z9QoIZT83I+5YNuv71KlTUlpaarq73fSXGjVqlBQUFJh9fdTubneQVnp+eHi4aYHXp7Ky0lykuhuAwKAMKNA8lg3UGqSVtqDr0n33MX3s2bOn1/F27dpJbGys5xxfcnNzTdB3b0lJSX75HQDUXwY0a0KyaUErfdR9d3lQAA7P+s7JyZHs7GzPvraoCdZA4ESEh8mytH6SNTHZjFnrqlh0dwM2a1HHx8ebx7KyMq/ndd99TB/Ly8u9jl+9etVkgrvP8SU6OtqMCdTdAAS+YIm7DChBGrBhoO7bt68Jtvn5+V4tXx17TklJMfv6WFFRIYWFhZ5zdu/eLbW1tWYsG4D/CpY8+3GRSQobt36PedR9fR5ACHV963znEydOeCWQHT161Iwx9+7dW5YsWSJPP/203HTTTSZwr1y50mRyuzPDBwwYIFOmTJHMzEwzhau6ulqysrJMRnhTM74BNB8FSwCHBOr//Oc/MmHCBM++e9w4PT3dTMF67LHHzFxrnRetLeexY8ea6Vft27f3fM+bb75pgvOkSZNMtvfs2bPN3GsAwSlYouPOdGUDITiPOpiYRw00nY5Ja3d3ff792AQz7gwgxOdRA7AmzdDuUE9BEp1mpccBtB0CNYBmoWAJEFiOnEcNoHXchUl0TPpyVY1pST88pi8FSwA/YIyaMWqgVYllFCwBmo8xagBNQsESwPro+gYcSAuT6FxonWalq1hpcljG2F+7rrW8JwDrIFADDkTBEsA+yPoGHKaxgiUtrdsNwD8I1IDDaPKXdnf7ohncehyAdRCoAYehYAlgLwRqwGEoWALYC8lkgANRsASwDwqeUPAEDkbBEiA4KHgCOAQFS4DQR9c3YEMULAGcg0AN2BAFSwDnIOsbsBkKlgDOQqAGbIaCJYCzEKgBm6FgCeAsBGrAZihYAjgLyWSADecyU7AEcA4KnlDwBDaeYkXBEiD0C57QogZsPMVKW+BJsTF+e48Ago8xaiAImGIFICQC9Zo1ayQsLMxr69+/v+f4lStXZNGiRdKtWzfp1KmTzJ49W8rKyoL6noGmYIoVgJAI1GrQoEFSUlLi2Q4cOOA5tnTpUvnggw9k+/btsm/fPjl79qzMmjUrqO8XaAqmWAFoKsuPUbdr107i4+OveV4H4F977TXZunWrTJw40Ty3efNmGTBggBw6dEhGjx4dhHcLNG+KVd0xareHx/RtcvY3gNBn+Rb1d999J4mJiXLDDTfI3Llz5fTp0+b5wsJCqa6ultTUVM+52i3eu3dvKSgoaPA1KysrTcZd3Q0I9ApWmt2dNSFZYqJ+Dcr6qPvuqVcAYPkW9ahRo2TLli3Sr18/0+29du1aGTdunBw7dkxKS0slKipKunbt6vU9cXFx5lhDcnNzzWsBwZxepedodnfWxORmz6MG4ByWDtRTp071fD106FATuPv06SPvvPOOdOjQocWvm5OTI9nZ2Z59bVEnJSW1+v3CWdpqBSumWAGwddd3Xdp6vvnmm+XEiRNm3LqqqkoqKiq8ztGsb19j2nVFR0ebCeZ1N6A5mF4FIFBsFagvXrwoJ0+elISEBBk+fLhERkZKfn6+53hRUZEZw05JSQnq+0ToY3oVgECxdNf3smXLZMaMGaa7W6derV69WiIiIuS+++4zpdcyMjJMF3ZsbKxpFS9evNgEaTK+EajpVdrd/XuaFKbHASDkA/WZM2dMUP7xxx+lR48eMnbsWDP1Sr9Wzz//vISHh5tCJ5rJnZaWJq+88kqw3zZspiX1spleBSBQWJSDRTkcq7VZ2+7vf/3gKblcVWNa0joHujmLagBwpp+bsSgHgZpA7VjPflzks+CIzmVuTtY2K1gB8GegtlUyGWDFrG339CrmQAPwBwI1HImsbQB2QaCGI7EoBgC7IFDDkfW23VnbvrAoBgArsfT0LMCfmdvuxS98ZW0DgFWQ9U3Wtzg9c5usbQCBRtY3HKGtMrfJ2gZgZYxRw7bI3AbgBARq2DIZTJG5DcAJSCaDbZPBqLcNwAkI1AgqDdJ1k8E0WLv3m5IMRuY2gFBH1jdZ30Gj3dy3PplX71KRR1be0eSynGRuA7ATsr5hizHmtkwGI3MbQKii6xtBG2N2J4PV16LW4wDgdGR9o03GmN3B1j3GrM83hjKeANA4AjVa3H3dFgVHtOWtVcS0Ba30Ufcp4wkAv6LrGy3uvm7KGLOu09wQfX3N7s6amGzO1+5u1nUGgP9HoEaLp0i15RizOxkMAOCNrm+HZ163pvuaMWYA8D9a1A7PvG5t9zUFRwDAvwjUFtLSoh2tqe7V2u5rxpgBwL/o+rZIi1jXVdYqXePW7zGPuq/P+zvzuq26ryk4AgD+QaAO8hhxa+cit0V1L6ZIAYB1hUzX98aNG2XDhg1SWloqw4YNk5dffllGjhxp+THixlrEOm2poVZtW2Re030NANYVEi3qt99+W7Kzs2X16tVy5MgRE6jT0tKkvLw8ID8/mC3itsy8pvsaAKwnJAL1c889J5mZmfLQQw/JwIEDZdOmTRITEyOvv/66z/MrKyvNyiV1t5Zq7Rixu0XsS1NbxHRdA0Dosn3Xd1VVlRQWFkpOTo7nufDwcElNTZWCggKf35Obmytr165tk5/f2ulN7hZx3azt5raI6boGgNBl+xb1Dz/8IDU1NRIXF+f1vO7reLUvGtTPnz/v2YqLi1v8863UIqbrGgBCj+1b1C0RHR1ttrZAixgA4E+2D9Tdu3eXiIgIKSsr83pe9+Pj4wPyHtqqOhf1rgEAIReoo6KiZPjw4ZKfny8zZ840z9XW1pr9rKysgLwHxogBAP5i+0CtdGpWenq6jBgxwsydfuGFF+TSpUsmCzyQaBEDANpaSATqe++9V86dOyerVq0yCWS33HKL7Nq165oEMwAA7CbM5XI1XlA6xOk86uuuu85kgHfp0iXYbwcAEOJ+bkbcsf30LAAAQhmBGgAACwuJMerWcvf+t6aUKAAATeWON00ZfSZQi8iFCxfMxUhKSmryRQYAoC3ij45VN4Rkst/mXZ89e1Y6d+5sLpoGbC0rSmJZ8z4dct1ahmvHdQs07rngXzdtSWu8SUxMNOtTNIQW9W+LePTq1ctckLCwX9eP1n8EAnXzcd1ajmvHdQs07rngXrfGWtJuJJMBAGBhBGoAACyMQP07uqrW6tWr22x1LafgunHtuOfsg/+v9rpuJJMBAGBhtKgBALAwAjUAABZGoAYAwMII1AAAWBiBuo6NGzfK9ddfL+3bt5dRo0bJZ599Frx/GZtYs2aNKRJTd+vfv3+w35bl7N+/X2bMmGGqEOk1eu+9966pUqTrqSckJEiHDh0kNTVVvvvuu6C9Xztdu3nz5l1zD06ZMkWcLjc3V2677TZTcbFnz54yc+ZMKSoq8jrnypUrsmjRIunWrZt06tRJZs+eLWVlZeJkuU24brfffvs199yCBQv89p4I1L95++23JTs726TeHzlyRIYNGyZpaWlSXl7ut4sfKgYNGiQlJSWe7cCBA8F+S5Zz6dIlc0/ph0Ff1q9fLy+99JJs2rRJDh8+LB07djT3n/4hdbrGrp3SwFz3Hty2bZs43b59+0wQPnTokOTl5Ul1dbVMnjzZXE+3pUuXygcffCDbt28352sp5VmzZomT7WvCdVOZmZle95z+H/YbF4yRI0e6Fi1a5LkaNTU1rsTERFdubi5XqAGrV692DRs2jGvUDPrfbseOHZ792tpaV3x8vGvDhg2e5yoqKlzR0dGubdu2cW0buHYqPT3dddddd3GdGlFeXm6u3759+zz3WGRkpGv79u2ec7755htzTkFBAdeznuum/vSnP7keeeQRV6DQohaRqqoqKSwsNN2Ndet/635BQYH/PiWFCO2i1W7JG264QebOnSunT58O9luylVOnTklpaanX/ac1gHX4hfuvafbu3Wu6Kfv16ycLFy6UH3/80W//XnZ1/vx58xgbG2se9W+ethbr3nc6bNW7d2/uuwaum9ubb74p3bt3l8GDB0tOTo5cvnxZ/IVFOUTkhx9+kJqaGomLi/O6OLr/7bff+u3ihwINJlu2bDF/ILX7Z+3atTJu3Dg5duyYGeNB4zRIK1/3n/sYpMFub+2u7du3r5w8eVIef/xxmTp1qgk2ERERXLrfVghcsmSJjBkzxgQW930XFRUlXbt25b5rxnVT999/v/Tp08c0UL788ktZsWKFGcd+9913xR8I1GgV/YPoNnToUBO49QZ+5513JCMjg6sLv5szZ47n6yFDhpj78MYbbzSt7EmTJvEvIGLGXPXDM/kjbXPd5s+f73XPaRKo3mv6QVHvvbZG17eI6b7QT96/z3bU/fj4+Da/6KFMP53ffPPNcuLEiWC/Fdtw32Pcf21Dh2D0/zT34K+ysrJk586dsmfPHs9yvu77Tof9KioqvK4ff/cavm6+aANF+eueI1CLmO6f4cOHS35+vleXh+6npKT45cKHqosXL5pPlfoJE02jXbb6R7Pu/acL1Gv2N/df8505c8aMUTv9HtTcOw02O3bskN27d5v7rC79mxcZGel132n3reaYOPm+czVy3Xw5evSoefTXPUfX9290alZ6erqMGDFCRo4cKS+88IJJx3/ooYf8cuFDxbJly8wcV+3u1qkdOr1Neyfuu+++YL81y32AqftpWxPI9D+3Jqho8o6Ogz399NNy0003mT8MK1euNONfOofT6Rq6drppXoTO/9UPO/oh8bHHHpPk5GQzvc3p3bZbt26V999/3+SLuPMdNFFR5+rrow5P6d8+vY5dunSRxYsXmyA9evRocapFjVw3vcf0+LRp08z8cx2j1mlu48ePN8MufhGw/HIbePnll129e/d2RUVFmelahw4dCvZbsrx7773XlZCQYK7ZH//4R7N/4sSJYL8ty9mzZ4+Z4vH7TacWuadorVy50hUXF2emZU2aNMlVVFQU7Ldt+Wt3+fJl1+TJk109evQwU4369OnjyszMdJWWlrqcztc1023z5s2ec3755RfXX/7yF9cf/vAHV0xMjOvuu+92lZSUuJxMGrlup0+fdo0fP94VGxtr/q8mJye7li9f7jp//rzf3hPLXAIAYGGMUQMAYGEEagAALIxADQCAhRGoAQCwMAI1AAAWRqAGAMDCCNQAAFgYgRoAAAsjUAMAYGEEagAALIxADQCAhbF6FoBmu/32281KQe3bt5e///3vZqnYBQsWyJo1a7iaQBujRQ2gRd544w3p2LGjWTd7/fr18uSTT0peXh5XE2hjrJ4FoEUt6pqaGvn3v//teU7XcZ84caKsW7eOKwq0IVrUAFpEu77rSkhIkPLycq4m0MYI1ABaJDIy0ms/LCxMamtruZpAGyNQAwBgYQRqAAAsjEANAICFkfUNAICF0aIGAMDCCNQAAFgYgRoAAAsjUAMAYGEEagAALIxADQCAhRGoAQCwMAI1AAAWRqAGAMDCCNQAAFgYgRoAALGu/wM1QX++kK3ocAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = np.arange(1, m + 1)\n", "\n", "plt.figure(figsize=(5, 3.2))\n", "plt.scatter(n, s_n, s=18)\n", "plt.xlabel(\"n\"); plt.ylabel(\"s_n\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ff4d014e", "metadata": {}, "source": [ "### Vectorization\n", "\n", "Just like R, NumPy applies arithmetic element by element — no loop required." ] }, { "cell_type": "code", "execution_count": 40, "id": "f9c2862f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 1.4142 1.7321 2. 2.2361 2.4495 2.6458 2.8284 3. 3.1623]\n", "[ 1 4 9 16 25 36 49 64 81 100]\n" ] } ], "source": [ "x = np.arange(1, 11)\n", "y = np.arange(1, 11)\n", "\n", "print(np.sqrt(x))\n", "print(x * y)" ] }, { "cell_type": "markdown", "id": "71e68303", "metadata": {}, "source": [ "> **Watch out:** a plain Python **list** is not a vector.\n", "> `[1,2,3] * 2` repeats the list; `np.array([1,2,3]) * 2` doubles each element." ] }, { "cell_type": "code", "execution_count": 41, "id": "fd76c262", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 3, 1, 2, 3]\n", "[2 4 6]\n" ] } ], "source": [ "print([1, 2, 3] * 2)\n", "print(np.array([1, 2, 3]) * 2)" ] }, { "cell_type": "markdown", "id": "6a2452eb", "metadata": {}, "source": [ "---\n", "# 8. Exercise — the `grade_it` function" ] }, { "cell_type": "code", "execution_count": 42, "id": "e3142c3c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'F'" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def grade_it(grade):\n", " if grade >= 90:\n", " return \"A\"\n", " elif grade >= 80:\n", " return \"B\"\n", " elif grade >= 70:\n", " return \"C\"\n", " else:\n", " return \"F\"\n", "\n", "\n", "grade_it(30)" ] }, { "cell_type": "markdown", "id": "d9e3a87a", "metadata": {}, "source": [ "---\n", "# 9. Exercise — a `for` loop that uses `grade_it`\n", "\n", "R:\n", "```r\n", "letter_grade <- rep(\"NA\", length(scores))\n", "for (i in 1:length(scores)) {\n", " letter_grade[i] <- grade_it(scores[i])\n", "}\n", "```" ] }, { "cell_type": "code", "execution_count": 44, "id": "1df8a451-36cb-4bde-bfe4-6a9dbee8b544", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[55, 82, 91, 40, 67, 78, 100, 62, 35, 88]\n" ] } ], "source": [ "print(scores)" ] }, { "cell_type": "code", "execution_count": 43, "id": "8958ee8a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['F', 'B', 'A', 'F', 'F', 'C', 'A', 'F', 'F', 'B']" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "letter_grade = [\"NA\"] * len(scores) # R: rep(\"NA\", length(scores))\n", "\n", "for i in range(len(scores)):\n", " letter_grade[i] = grade_it(scores[i])\n", "\n", "letter_grade" ] }, { "cell_type": "code", "execution_count": 45, "id": "1168a280", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['F', 'B', 'A', 'F', 'F', 'C', 'A', 'F', 'F', 'B']" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# More Pythonic: append as you go, no placeholder needed\n", "letter_grade = []\n", "for score in scores:\n", " letter_grade.append(grade_it(score))\n", "\n", "letter_grade" ] }, { "cell_type": "markdown", "id": "06b51982", "metadata": {}, "source": [ "---\n", "# 10. `sapply()` → `map`, comprehensions, and `.apply()`\n", "\n", "R:\n", "```r\n", "x <- 1:10\n", "sapply(x, sqrt)\n", "```\n", "\n", "Python offers **four** natural replacements. Pick whichever reads best." ] }, { "cell_type": "code", "execution_count": 46, "id": "189d973d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[np.float64(1.0), np.float64(1.4142135623730951), np.float64(1.7320508075688772), np.float64(2.0), np.float64(2.23606797749979), np.float64(2.449489742783178), np.float64(2.6457513110645907), np.float64(2.8284271247461903), np.float64(3.0), np.float64(3.1622776601683795)]\n", "[np.float64(1.0), np.float64(1.4142135623730951), np.float64(1.7320508075688772), np.float64(2.0), np.float64(2.23606797749979), np.float64(2.449489742783178), np.float64(2.6457513110645907), np.float64(2.8284271247461903), np.float64(3.0), np.float64(3.1622776601683795)]\n", "[1. 1.4142 1.7321 2. 2.2361 2.4495 2.6458 2.8284 3. 3.1623]\n", "[1.0, 1.4142135623730951, 1.7320508075688772, 2.0, 2.23606797749979, 2.449489742783178, 2.6457513110645907, 2.8284271247461903, 3.0, 3.1622776601683795]\n" ] } ], "source": [ "x = list(range(1, 11))\n", "\n", "print([np.sqrt(i) for i in x]) # 1. list comprehension (most common)\n", "print(np.sqrt(x)) # 3. vectorized ufunc (fastest)\n", "print(pd.Series(x).apply(np.sqrt).tolist()) # 4. pandas .apply()" ] }, { "cell_type": "code", "execution_count": 47, "id": "3676e527-6b12-4bce-8b9c-0aff12bbcebb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[np.float64(1.0), np.float64(1.4142135623730951), np.float64(1.7320508075688772), np.float64(2.0), np.float64(2.23606797749979), np.float64(2.449489742783178), np.float64(2.6457513110645907), np.float64(2.8284271247461903), np.float64(3.0), np.float64(3.1622776601683795)]\n" ] } ], "source": [ "print(list(map(np.sqrt, x))) # 2. map()" ] }, { "cell_type": "markdown", "id": "3e9c5626", "metadata": {}, "source": [ "### Exercise — Build a grading machine with \"sapply\"" ] }, { "cell_type": "code", "execution_count": 48, "id": "e41ee3af", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['F', 'B', 'A', 'F', 'F', 'C', 'A', 'F', 'F', 'B']" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scores = [55, 82, 91, 40, 67, 78, 100, 62, 35, 88]\n", "\n", "letter_grade = [grade_it(s) for s in scores] # R: sapply(scores, grade_it)\n", "letter_grade" ] }, { "cell_type": "code", "execution_count": 49, "id": "919bf898", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['F', 'B', 'A', 'F', 'F', 'C', 'A', 'F', 'F', 'B']" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# the map() version\n", "list(map(grade_it, scores))" ] }, { "cell_type": "code", "execution_count": 50, "id": "451ac5e1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
scoreletterpass
055FFAIL
182BPASS
291APASS
340FFAIL
467FPASS
578CPASS
6100APASS
762FPASS
835FFAIL
988BPASS
\n", "
" ], "text/plain": [ " score letter pass\n", "0 55 F FAIL\n", "1 82 B PASS\n", "2 91 A PASS\n", "3 40 F FAIL\n", "4 67 F PASS\n", "5 78 C PASS\n", "6 100 A PASS\n", "7 62 F PASS\n", "8 35 F FAIL\n", "9 88 B PASS" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# the pandas version -- and a tidy gradebook to go with it\n", "gradebook = pd.DataFrame({\"score\": scores})\n", "gradebook[\"letter\"] = gradebook[\"score\"].apply(grade_it)\n", "gradebook[\"pass\"] = np.where(gradebook[\"score\"] >= 60, \"PASS\", \"FAIL\")\n", "gradebook" ] }, { "cell_type": "code", "execution_count": 52, "id": "a276f084-b465-40a0-bbe7-0f521ad29c41", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "letter\n", "F 5\n", "B 2\n", "A 2\n", "C 1\n", "Name: count, dtype: int64" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gradebook[\"letter\"].value_counts()" ] }, { "cell_type": "code", "execution_count": 51, "id": "8649fe0b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "letter\n", "A 2\n", "B 2\n", "C 1\n", "F 5\n", "Name: count, dtype: int64" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gradebook[\"letter\"].value_counts().sort_index() ### sort_index: sort by index alphabetically," ] }, { "cell_type": "markdown", "id": "14419ae6", "metadata": {}, "source": [ "---\n", "# 11. `apply()` over a matrix → `axis=`\n", "\n", "R:\n", "```r\n", "mat <- matrix(1:20, ncol = 2)\n", "apply(mat, 2, mean) # column means\n", "apply(mat, 1, mean) # row means\n", "```\n", "\n", "> **R vs. Python** R fills a matrix **column by column**; NumPy's `reshape` fills **row by row**.\n", "> Use `order=\"F\"` to reproduce R's layout." ] }, { "cell_type": "code", "execution_count": 53, "id": "94c325b9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1, 11],\n", " [ 2, 12],\n", " [ 3, 13],\n", " [ 4, 14],\n", " [ 5, 15],\n", " [ 6, 16],\n", " [ 7, 17],\n", " [ 8, 18],\n", " [ 9, 19],\n", " [10, 20]])" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat = np.arange(1, 21).reshape(10, 2, order=\"F\") # same as R's matrix(1:20, ncol=2)\n", "mat" ] }, { "cell_type": "code", "execution_count": 54, "id": "1e4031c9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cMean: [ 5.5 15.5]\n", "rMean: [ 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.]\n" ] } ], "source": [ "cMean = mat.mean(axis=0) # R: apply(mat, 2, mean) -- collapse the rows -> column means\n", "rMean = mat.mean(axis=1) # R: apply(mat, 1, mean) -- collapse the cols -> row means\n", "\n", "print(\"cMean:\", cMean)\n", "print(\"rMean:\", rMean)" ] }, { "cell_type": "markdown", "id": "bcfcbbc1", "metadata": {}, "source": [ "> **Remember the axis rule:** `axis=0` works *down* the rows (giving one number per column),\n", "> `axis=1` works *across* the columns (one number per row).\n", "> That is the mirror image of R's `apply(mat, 1, ...)` / `apply(mat, 2, ...)` numbering." ] }, { "cell_type": "code", "execution_count": 55, "id": "53651ec4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 5.5 15.5]\n", "[ 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.]\n" ] } ], "source": [ "# np.apply_along_axis() is the literal translation of apply(), and works\n", "# for any function -- not just the ones NumPy has built in.\n", "print(np.apply_along_axis(np.mean, 0, mat))\n", "print(np.apply_along_axis(np.mean, 1, mat))" ] }, { "cell_type": "code", "execution_count": 56, "id": "bbf96a28", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[5.5, 15.5]\n", "[6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0]\n" ] } ], "source": [ "# pandas keeps R's 1/2 numbering readable with axis names\n", "df = pd.DataFrame(mat, columns=[\"V1\", \"V2\"])\n", "print(df.apply(np.mean, axis=\"index\").tolist()) # column means\n", "print(df.apply(np.mean, axis=\"columns\").tolist()) # row means" ] }, { "cell_type": "markdown", "id": "293d94f7", "metadata": {}, "source": [ "---\n", "# 12. Putting it together — a simulation\n", "\n", "**Task**\n", "1. Generate 30 random numbers from a normal distribution.\n", "2. Randomly assign 15 to control and 15 to treatment.\n", "3. Run a *t*-test and report the p-value." ] }, { "cell_type": "code", "execution_count": 57, "id": "590f3da1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.9891, -0.3678, 1.2879, 0.194 , 0.9202])" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rng = np.random.default_rng(123) # R: set.seed(123)\n", "\n", "obs = rng.normal(loc=0, scale=1, size=30) # R: rnorm(30, mean=0, sd=1)\n", "obs[:5]" ] }, { "cell_type": "code", "execution_count": 58, "id": "9c68b119", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['B'], dtype=' **R vs. Python** R's negative index `obs[-icontrol]` means \"everything *except*\".\n", "> In Python a negative index counts from the end, so use a boolean mask or `np.setdiff1d()` instead." ] }, { "cell_type": "code", "execution_count": 60, "id": "8f287383", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TtestResult(statistic=np.float64(-0.8220123018732333), pvalue=np.float64(0.41801502605208213), df=np.float64(28.0))\n", "p-value: 0.41801502605208213\n" ] } ], "source": [ "icontrol = rng.choice(30, size=15, replace=False) # R: sample(1:30, 15)\n", "\n", "mask = np.zeros(30, dtype=bool)\n", "mask[icontrol] = True\n", "\n", "control = obs[mask] # R: obs[icontrol]\n", "trt = obs[~mask] # R: obs[-icontrol] (~ means \"not\")\n", "\n", "res = stats.ttest_ind(control, trt)\n", "print(res) # R: str(res)\n", "print(\"p-value:\", res.pvalue) # R: res$p.value" ] }, { "cell_type": "markdown", "id": "6c1c43b1", "metadata": {}, "source": [ "### Do it 1000 times — with a `for` loop\n", "\n", "R:\n", "```r\n", "sim <- 1000; n <- 15\n", "pvalue <- numeric(sim)\n", "for (i in 1:sim) { ... pvalue[i] <- res$p.value }\n", "```" ] }, { "cell_type": "code", "execution_count": 61, "id": "761297a9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.2491, 0.0894, 0.9008, 0.2565, 0.0359])" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sim = 1000\n", "n = 15\n", "rng = np.random.default_rng(1)\n", "\n", "pvalue = np.zeros(sim) # R: numeric(sim)\n", "\n", "for i in range(sim):\n", " obs = rng.normal(size=2 * n)\n", " icontrol = rng.choice(2 * n, size=n, replace=False)\n", " mask = np.zeros(2 * n, dtype=bool)\n", " mask[icontrol] = True\n", " res = stats.ttest_ind(obs[mask], obs[~mask])\n", " pvalue[i] = res.pvalue\n", "\n", "pvalue[:5]" ] }, { "cell_type": "markdown", "id": "eb9107a0", "metadata": {}, "source": [ "### The same thing as a function\n", "\n", "Note that `iter` is never used inside the function — it is just a placeholder so that\n", "`map`/`sapply` has something to iterate over. This is a common pattern in simulations.\n", "\n", "> **Naming note:** `iter` is a built-in Python function, so we call the argument `i` instead." ] }, { "cell_type": "code", "execution_count": 62, "id": "f3a01f96", "metadata": {}, "outputs": [], "source": [ "def getpvalue(i, n, rng):\n", " obs = rng.normal(size=2 * n)\n", " icontrol = rng.choice(2 * n, size=n, replace=False)\n", " mask = np.zeros(2 * n, dtype=bool)\n", " mask[icontrol] = True\n", " res = stats.ttest_ind(obs[mask], obs[~mask])\n", " return res.pvalue" ] }, { "cell_type": "code", "execution_count": 64, "id": "78c7b3cf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n" ] } ], "source": [ "# R: pvalue_sapply <- sapply(1:1000, getpvalue, n = 15)\n", "rng = np.random.default_rng(1)\n", "pvalue_listcomp = np.array([getpvalue(i, n=15, rng=rng) for i in range(1000)])\n", "\n", "# R's replicate() equivalent\n", "rng = np.random.default_rng(1)\n", "pvalue_map = np.array(list(map(lambda i: getpvalue(i, 15, rng), range(1000))))\n", "\n", "# R: identical(pvalue_sapply, pvalue_replicate)\n", "print(np.array_equal(pvalue_listcomp, pvalue_map))" ] }, { "cell_type": "code", "execution_count": 65, "id": "a9a8a159", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean p-value: 0.5038851300705502\n", "proportion < 0.05: 0.055 (should be about 0.05)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE1CAYAAAAoBaPAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMqdJREFUeJzt3QmcTfX/+PH3jGHGOvatrCFSsoWxJKKJEqWyZcuXFhQS+VJaZCtbZSmVpUIURWoUQhiFIvHNniWGFmPNEOf/eH/+j3N/9965w9yZa+6Zua/n43GMe+65537u55573ueznjDLsiwBAACOFB7sBAAAgJQRqAEAcDACNQAADkagBgDAwQjUAAA4GIEaAAAHI1ADAOBgBGoAAByMQA0AgIMRqBEUYWFh8uKLL7oez5w506z77bff0r3vO+64wyzBEMjPYdN96T5130ibVatWmTzUv8gYvn6Hx44dkwcffFAKFSpkvo+JEyfydaQCgToDnDlzRoYPHy533323FCxY8Kon3f/9739m2zx58pjtO3fuLH/88Uey7S5fvixjx46VcuXKSVRUlFSrVk3mzp2brn0ieObMmcOJC1la//79ZdmyZTJkyBD54IMPzDkJVxeRim2QTn/++ae8/PLLUrp0abn11luveFV/+PBhuf322yU6OlpGjhxpgvzrr78u27Ztkx9++EFy5Mjh2nbo0KEyevRo6dmzp9x2223y+eefS8eOHc2FQPv27dO0z6zg66+/lswaqH/55Rfp16+fx/oyZcrIP//8I9mzZw9a2oBA/A5XrlwprVu3loEDB5KhfiBQZ4ASJUrI0aNHpXjx4rJp0yYTVFOigfTs2bOyefNmE9hVnTp1pHnz5qYU3qtXL7Pu999/l3Hjxknv3r3lrbfeMuv+85//SOPGjeXZZ5+Vhx56SLJly+bXPrOKrHbhoRdeWmMCeNN7Kp0/f15y5syZKX6Hx48fl/z58wfsPc6fP2/eJzw8a1cOZ+1P5xCRkZEmSKfGp59+Kvfee68roKpmzZpJpUqVZP78+a51Wnq+ePGiPPnkkx4n9CeeeMKUoOPj4/3e55XaR7UE/s4778gNN9xgPo9ebGzcuDFVbcPdunWTsmXLSiAkJCRI9+7d5frrrzfp0IsgvUJ3bxP2TofdPqmf9aWXXpLrrrtO8ubNa9rKTp48KUlJSaYUW7RoUdM0oPvXdd554Ku5wrut3Rf9ru655x4pWbKkSbPm4SuvvCKXLl3ySPPSpUvlwIEDZp+62HmW0vtr6aRRo0aSO3duc/LTfNAmDneaNn3tnj17zPeg22nNin7Gc+fOeWz7zTffSMOGDc02mg833nij/Pe///XY5uDBg/Lrr79e8fO65/nHH39s9qHHv6bzvvvuk0OHDl3xtXoxq6+dNWtWsue02lSf++KLL8xjzS/9DWhaNVhp26depKamj4Dmr+aJN1/HsR4P2nxVoUIF8x2WKlVKBg0a5HGcpDYPfZkxY4Y0bdrUHIO6/5tuukmmTp3qM836W9Z8qF27tvnMb7/9tnkuMTHRHMeaNt2HpnXMmDGmiexqUjqOvfPI7oOxbt06GTBggBQpUsR8r/fff3+ypjT3fLRfpxcWkydPdh3jtn379pnvTZvlcuXKJfXq1TO/B1/H1Lx582TYsGHmd6zbnjp1yqRR81uPT80f/b8+r++ltPZQ81fTqjVUWnuVmVCidhAtJesVp/4AvWkJ+Msvv3Q9/umnn8xBV6VKlWTb2c/rCcOffV6JHtinT5+Wxx57zPxYtG38gQceMD+wjKySbdu2rWzfvl369u1rTiL62fTkqD/Qq10MjBo1ypzYnnvuORO43nzzTZN2vRo/ceKEOVFt2LDBnFS03f+FF14ISJp1f3ri0BOb/tUAq/vWE8xrr73masbQiwa9yJowYYJZp9umZPny5dKiRQspX768SbdWjevnadCggfz444/J8uLhhx82n0nzQJ9/9913TVDQE7nSPNUTnPZz0GYaPdFrHukJ2V2XLl1k9erV5oSbGq+++qo5XgYPHmy+K+08pBeJW7ZsSbEUqMeqfi69sOratavHcxr4CxQoILGxseaxXiyuX7/eNPXoxZsGaA1wGiB27NhhTuTppYFOLzDWrl1rap/0N6cnfv2edu3aJZ999plfeeiLprlq1armfSIiImTJkiXmAkTfW2vN3O3cuVM6dOhgfova7KUXA3rRpbVp+nvX9XpRrvmibcFamxfoTlv6+9PvQS9eNM91/3369DHfjy/a9KZt0to3Rmvy9Dhy72BWv3598xmeeuopc7GlF2maF5988om5CHCnF7laitbqc71QskvueuGrvwl9Lz0/ffTRRyZNep7U31enTp3MOWvatGnm/WNiYsxvIlPQ+1Ej42zcuFHPcNaMGTNSfG727NnJnnv22WfNc+fPnzeP77nnHqt8+fLJtjt79qzZ7rnnnvN7n77s37/fbFOoUCHr77//dq3//PPPzfolS5a41jVu3Ngs3rp27WqVKVPGY52+dvjw4a7Hmh+6Tt8vJSdOnDDbvPbaaylu4ysd3377rXndzTffbF24cMG1vkOHDlZYWJjVokULj9fHxMR4pNfOA1/fWWo+x7lz55K97rHHHrNy5crlkff6nXrnU0rvX716dato0aLWX3/95Vq3detWKzw83OrSpYtrnaZNX/voo4967PP+++8336ltwoQJZrs//vjDuhLN19ScNuw8v+6666xTp0651s+fP9+snzRp0hVfP2TIECt79uwex1xSUpKVP39+j8/iK2/j4+OTHfN2evSvTfNaj82rHT8ffPCBydfvvvvOY7tp06aZfa5bt86vPPTF1+eIjY1N9hvXNOt7xMXFeax/5ZVXrNy5c1u7du3yWK/ngWzZslkHDx684vt7H8cp5ZF9fDdr1sy6fPmya33//v3N+yQmJl7xfKCv7d27t8e6fv36mfXu+Xv69GmrXLlyVtmyZa1Lly55fIeaJ975pWnU50aOHOlxvsiZM6f5jc+bN8+1/tdff03x8zoVVd8OoqUipVfi3uw2Snsb/Zva7VK7zytp166duYK2aZWr0hJ1RtESmF49axWYloD9pVfR7qX/unXrmpLho48+6rGdrtfq2X///Tdg6bZprYR2LtT80xJEaqqRvWkJSUukWt2nVYU2LclpacVXLcnjjz/u8Vjf/6+//jKlemW3G2o1/ZWqSjXvU1uatvNcmxls2tygzRVXq8nR402bdhYuXOjROUmrd/U5X3mr2+tn0ipf/TxacxAICxYsMKXoypUrm+/OXrQqVX377bd+5aEv7p9Da1Z0/1pC1t+XPnanpUC7RsE9jfqd6m/UPY1ae6ElzTVr1kggac2Ce9W1vre+jzZF+EuPBa3d0xpAm9Ym6XtoaX3Hjh0e22stS0q1MdpPx6bfh9Y2aIlaa5Rsuk6fy8hzV3oRqB3EPvi8273sThPu2+jf1G6X2n1eiXv7trKDdloCZlrpxYZW1X711VdSrFgxVxWXtlunhvdn0LZapW163uv1ROt9gkwrrRLV6jvdb758+Uy73iOPPGKeS8t72CdDPeF404CiJ2jtPOjP96fBT6vN9USneatVyVr17G/A8VaxYkWPx3py10BqtyHrCAT9/uzFbufU0REaGN2rUvX/hQsXdgVI+yJTmxHsdll9XvNXA3qgvr/du3eb71D3675oHw+lVfrpzUOtHtegavc30P3bbdu+ArWvNMbFxSVLo+7TPY2BEsjzgR7PKR3L9vPuUqqu1oKHfmZ3+pvTJhH3iwp7fUaeu9KLNmoH0ZKGXWLypuu09GSXjHVbvZLX0o37QWi/Vjsu+bvPK7F7kHtzL13ZnUW8uXeaSi/tLNOqVSvTLqgdap5//nnT7qrtvjVq1EjTZ7jaZ/P+kfvzuTRgaMlIA7S2W2pHMj2haGlP223TGwhT62qfUS/WtNSlx5R24tGTvgZGDYpakk3p9emlnRS1g59NO/rYQVwDn7Zx64WHlsoXL15s2ma1Dde9rVQ7YulxoW2OegK2hydeLW+v9L26f17dzy233CLjx4/3ub19oZfWPNy7d6/ceeed5sJE30P3pzVHWtLUdnDvz+Hrwlq30doU7eDmi31R4a+UjvHUnA+ulZwpFCzS+vvODAjUDqK9FPWKUHu9etPxztWrV3c91v9rhyDt5as9RG3ff/+963l/95leelXtqzopLdVhV6LB7plnnjGLliT0M+hQtQ8//FCuBbu0oEHX38+lVcVaHatVuFoDYNu/f3+qA4c3DWZ2pyJvWpWupUotmflLO9VpwNBFA4YO69NOOBp47JKZv/T78T45agcrraa3q8bdqzzdT8IaqDWI66gFLaFqNb37/ABKOxtpVah+/+41Rd7fVUrfq6/t9HvVzmzux9vWrVtNvlztO0pLHmrHMa3x0gsR95KqXaWeGppGrZ1I6/fkKy8uXLjg8wI/0PR4TulYtp8PdVR9O4z2atahJ+5DWFasWGF6l+rwBZsOxdH21ilTpnicBLVHowZn7UXp7z7TS08W+uNyH6ahJ7jU9HpNDW3Ttavr3d9TS1u+qvYDRUvDGvy82/nc8z4l9tW8+9W7ngB9vVaDa2qqa7WWRC9OtGes+8lVJ0vRklvLli3FX3///XeydfZFnHvepnZ4lm327NmmXd49sOrJX3vnKg2IGlzsRauO3as+tSSrpVJd9HO7X+zY+etdMtLe76mp7dBjR3v56/dh8/6dKG3f1N7U06dPT7YPrXq3mxlSm4epOUb0ONCagtTSNOqQTK1l8qbHyNX6W2heeB/fOhwzkLVhKdHjVQsN7kNKNU/1/XX0wk1uBZFQRYk6g+ikJPqDOXLkiOsqWofi2NV3dnuptktpx5AmTZrI008/ba6SdQiPnrB07KtN2120uk+f0040Oq5Zq4O/++47MyzBvbontftML+2UpaUI7ejSo0cP0y6mFw467MTutJQeemGhJRU9KemPV6tAFy1aZIZ3eJe0Ak3bHXUWOP2rw4f0pKbpuRq9YNLSipb6dOiJlsh0mIqvardatWqZgKTDuPT71A41Ws3vi35/Guy0ulfz2h6epcfR1cZ1+6LV8vqZdLy3lmD0u9OLCT3O3Eu8/g7P0qYVfb0eZ/o96TAebaPWYUWpoaVqbYPW5gL9nN4TW+hwKM1P/dx6TOjJXoeu6RCfq9HvUi8cdBpLPaa0ClprZTRoudMhRdrWrB3ytJSrFxMawPSCRdfbY5pTm4fe7rrrLlPVrd+1Dq3S36deFOjwudSWaHWSIy2Ra35oJ0M9ljTY6TAy/YzanKAXm1fKC/18elGvVeh6ga2f60qvCRQdLqlTH+vxrL8RPWb0IlRrnbQ2JTyLT2aSKsHudh4q7GEVvhbvIUm//PKLddddd5nhOzocpVOnTlZCQkKyfeqwBR2OoPvOkSOHVbVqVevDDz/0+f6p3WdKQ4N8DYnyNcRB31+HT2h6dAjRsmXLAjY8688//zRDOypXrmyGokRHR1t169Y1Q35SMzxrwYIFHtvZ76lD2NzZQ5rch9nocJAePXqY98ybN6/18MMPW8ePH0/V59DhO/Xq1TNDRUqWLGkNGjTI5Iv3cKEzZ85YHTt2NN+PPmfnWUrDw5YvX241aNDA7DdfvnxWq1atrB07dlz1s/hK54oVK6zWrVub9Ol3p391+Jr3cB9/h2fNnTvXDLXSoWSaTh2CduDAASu1du/e7fqdrF27NtnzOgSne/fuVuHCha08efKYIU06/MZ7WJGv4Vlq3LhxZghZZGSkyctNmzb5HFakw/rGjBljfmO6bYECBaxatWpZL730knXy5Em/8tCXxYsXW9WqVbOioqLMkCR9r/fffz/ZsaSfS/PQFx3SpHldoUIF8/6aJ/Xr17def/11j2GJvui5ZPDgweY1eo7QfNyzZ0+Kw7O8fzO+8je1w7PU3r17rQcffNAc+5oHderUsb744guf77HA63esNI16TvCm76/fmbcr5aMThek/wb5YAJC1aNu81uBoTY4OyQKQdtQpAADgYARqAAAcjEANAICD0UYNAICDUaIGAMDBsvw4ap1aT8cu66QYqZ35CQCAa0kHXOlkQDrd89XGimf5QK1B2vumCwAAOIHOhKeT4oR0oLZvsaeZoVNBAgAQbDpboxYi3W8DG7KB2q7u1iBNoAYAOElqmmTpTAYAgIMRqAEAcDACNQAADkagBgDAwQjUAAA4GIEaAAAHI1ADAOBgBGoAABws6IH6999/l0ceeUQKFSokOXPmlFtuuUU2bdrkMR/qCy+8ICVKlDDPN2vWTHbv3h3UNAPI2i5dthy5L4SmoM5MduLECWnQoIE0adJEvvrqKylSpIgJwgUKFHBtM3bsWHnjjTdk1qxZUq5cOXn++eclNjZWduzYIVFRUcFMPoAsKlt4mDw97yfZc/xMuvZToWgemdS+RsDShdAU1EA9ZswYM9fpjBkzXOs0GLuXpidOnCjDhg2T1q1bm3WzZ8+WYsWKyWeffSbt27cPSroBZH0apLcfORXsZADBrfpevHix1K5dWx566CEpWrSo1KhRQ6ZPn+56fv/+/ZKQkGCqu23R0dFSt25diY+P97nPpKQkM9m5+wIAQGYV1EC9b98+mTp1qlSsWFGWLVsmTzzxhDz11FOmmltpkFZagnanj+3nvI0aNcoEc3vhFpcAgMwsqIH68uXLUrNmTRk5cqQpTffq1Ut69uwp06ZNS/M+hwwZIidPnnQtentLAAAyq6AGau3JfdNNN3msq1Klihw8eND8v3jx4ubvsWPHPLbRx/Zz3iIjI123tOTWlgCAzC6ogVp7fO/cudNj3a5du6RMmTKujmUakFesWOF6Xtucv//+e4mJicnw9AIAkNGC2uu7f//+Ur9+fVP1/fDDD8sPP/wg77zzjlnsG2r369dPRowYYdqx7eFZJUuWlDZt2gQz6QAAZP1Afdttt8miRYtMu/LLL79sArEOx+rUqZNrm0GDBsnZs2dN+3ViYqI0bNhQ4uLiGEMNAAgJYZYOVs7CtKpce39rxzJtswaA1Ljnje/SPY66asl8svSpRmQ40hWbgj6FKAAASBmBGgAAByNQAwDgYARqAAAcjEANAICDEagBAHAwAjUAAA5GoA6iS5ctR+8PuNYCecxy/AcX32UWnZks1GULD5On5/1kblCfXhWK5pFJ7WsEJF1AZvsNcPwHH9/ltUOgDjI9QaV39iMgM+M3kHXwXV4bVH37ieo1AKlVJE8kTVxIN0rUQareuePGIvJsbOV07QOAs+XLGUETF9KNQB2k6p0biuRO1+sBZB5UCSM9qPoGAMDBCNQAAMegXT85qr4BAI5Bu35yBGoAgOPQrv9/qPpGpsJsbsio6lLAKShRI1NhNjdc6+pShk7CaQjUyHSoEsO1PDYYOgmnoeobAAAHI1DDJ+6EAwDOQNU3fOJOOADgDARqpIi2YAAI8arvF198UcLCwjyWypX/70YV58+fl969e0uhQoUkT5480rZtWzl27FgwkwwAQGi1UVetWlWOHj3qWtauXet6rn///rJkyRJZsGCBrF69Wo4cOSIPPPBAUNMLAEBIVX1HRERI8eLFk60/efKkvPfeezJnzhxp2rSpWTdjxgypUqWKbNiwQerVqxeE1AIAEGIl6t27d0vJkiWlfPny0qlTJzl48KBZv3nzZrl48aI0a9bMta1Wi5cuXVri4+NT3F9SUpKcOnXKYwEAILMKaqCuW7euzJw5U+Li4mTq1Kmyf/9+adSokZw+fVoSEhIkR44ckj9/fo/XFCtWzDyXklGjRkl0dLRrKVWqVAZ8EgBAKEwveykIU9UGteq7RYsWrv9Xq1bNBO4yZcrI/PnzJWfOnGna55AhQ2TAgAGux1qiJlgDQGjKF8DpZSsUzSOT2teQkGujdqel50qVKsmePXukefPmcuHCBUlMTPQoVWuvb19t2rbIyEizAACQFYabBr2N2t2ZM2dk7969UqJECalVq5Zkz55dVqxY4Xp+586dpg07JiYmqOl0Iu4eBABZU1BL1AMHDpRWrVqZ6m4dejV8+HDJli2bdOjQwbQv9+jRw1RjFyxYUPLlyyd9+/Y1QZoe38lx9yAAyJqCGqgPHz5sgvJff/0lRYoUkYYNG5qhV/p/NWHCBAkPDzcTnWhv7tjYWJkyZUowk+x43D0IALKWoAbqefPmXfH5qKgomTx5slkAAAhFjmqjBhB4mX1oCq4dvs/MwVG9vgEEXmYfmgLnHxt33FhEno39v/s0ILAI1EAIyMxDU3Bt0a/F+aj6RsjKCjMWAcj6KFEjZGWFGYsAZH0EaoQ8qoUBOBlV3wAAOBiBGgAAByNQAwDgYARqAAAcjEANAICDEagBAHAwAjUAAA5GoAYAwMEI1IADMR0pABszkwEOxF2NANgI1IBDcVcjAIqqbwAAHIxADQCAgxGoAQBwMAI1AAAORqAGAMDBCNQAAGSlQL1v375rkxIAAJD+QF2hQgVp0qSJfPjhh3L+/HkJlNGjR0tYWJj069fPtU7337t3bylUqJDkyZNH2rZtK8eOHQvYewKBUiRPZJafTSwUPqPT8R2EJr8nPPnxxx9lxowZMmDAAOnTp4+0a9dOevToIXXq1ElzIjZu3Chvv/22VKtWzWN9//79ZenSpbJgwQKJjo427/fAAw/IunXr0vxewLWQL2dElp9NLJCf0cmf08lC4ThDAAJ19erVZdKkSTJu3DhZvHixzJw5Uxo2bCiVKlWSRx99VDp37ixFihRJ9f7OnDkjnTp1kunTp8uIESNc60+ePCnvvfeezJkzR5o2bWrW6QVClSpVZMOGDVKvXj1/kw5cc6Ewm1ggPmNm+JxOFgrHGQLQmSwiIsKUbrW0O2bMGNmzZ48MHDhQSpUqJV26dJGjR4+maj9atX3PPfdIs2bNPNZv3rxZLl686LG+cuXKUrp0aYmPj09xf0lJSXLq1CmPBQCAkAvUmzZtkieffFJKlCgh48ePN0F679698s0338iRI0ekdevWV93HvHnzTFX6qFGjkj2XkJAgOXLkkPz583usL1asmHkuJbovrSa3F71wQPDQpgYAGVz1rUFZq6B37twpLVu2lNmzZ5u/4eH/P+aXK1fOVIeXLVv2ivs5dOiQPP300yawR0VFSaAMGTLEtJ/btERNsA4e2tQAIIMD9dSpU01bdLdu3Uxp2peiRYua9uUr0art48ePS82aNV3rLl26JGvWrJG33npLli1bJhcuXJDExESPUrX2+i5evHiK+42MjDQLnIU2NQDIoEC9e/fuq26jVdZdu3a94jZ33nmnbNu2zWNd9+7dTTv04MGDTSk4e/bssmLFCjMsS2kp/uDBgxITE+NvsgEACI1ArdXeOqb5oYce8livncrOnTt31QBty5s3r9x8880e63Lnzm3GTNvrddiXVmMXLFhQ8uXLJ3379jVBmh7fAIBQ4XdnMu2sVbhwYZ/V3SNHjpRAmjBhgtx7772mRH377bebKu+FCxcG9D0AAMhSJWqtetYOY97KlCljnkuPVatWeTzWTmaTJ082CwAAocjvErWWnH/++edk67du3WqqrQEAQBADdYcOHeSpp56Sb7/91vTS1mXlypVmqFX79u0DmDQAAOB31fcrr7wiv/32m+m1rbOTqcuXL5vZyALdRg0AQKjzO1Dr0KuPP/7YBGyt7s6ZM6fccsstpo0aAAAEOVDb9CYcugAAAAcFam2T1ilCdSISnVlMq73daXs1AAAIUqDWTmMaqPWOVzoxSVhYWICSAgAA0h2o9Y5X8+fPNzfiAAAADhuepZ3JKlSocG1SAwAA0heon3nmGZk0aZJYluXvSwEAwLWu+l67dq2Z7OSrr76SqlWrmjtcuWMubgAAghio9d7Q999/fwCTAAAAAnqbSwAA4NA2avXvv//K8uXL5e2335bTp0+bdUeOHJEzZ84EOn0AAIQ0v0vUBw4ckLvvvtvc0jIpKUmaN28uefPmlTFjxpjH06ZNuzYpBQAgBIWnZcKT2rVry4kTJ8w83zZtt9bZygAAQBBL1N99952sX7/ejKd2V7ZsWfn9998DmDQAAOB3iVrn9tb5vr0dPnzYVIEDAIAgBuq77rpLJk6c6Hqsc31rJ7Lhw4czrSgAAMGu+h43bpzExsbKTTfdJOfPn5eOHTvK7t27pXDhwjJ37txApw8AgJDmd6C+/vrrZevWrebmHD///LMpTffo0UM6derk0bkMAAAEIVCbF0VEyCOPPBKAtwcAAAEN1LNnz77i8126dPF3lwAAIFCBWsdRu7t48aKcO3fODNfKlSsXgRoAgGD2+taJTtwXbaPeuXOnNGzY0O/OZFOnTpVq1apJvnz5zBITE2PuymXTzmq9e/eWQoUKSZ48eaRt27Zy7Ngxf5MMAEBozfXtrWLFijJ69Ohkpe3UdEzT123evFk2bdokTZs2ldatW8v27dvN8/3795clS5bIggULZPXq1WY+8QceeCAQSQYAIOt2JvO5o4gIE0j90apVK4/Hr776qillb9iwwQTx9957T+bMmWMCuH3nripVqpjn69WrF6ikAwCQdQL14sWLPR5bliVHjx6Vt956Sxo0aJDmhOhsZ1pyPnv2rKkC11K2tn83a9bMtU3lypWldOnSEh8fn2Kg1huD6GI7depUmtMEAECmC9Rt2rTxeKwzkxUpUsSUenUyFH9t27bNBGZtj9Z26EWLFpnJVLZs2WI6qOXPn99j+2LFiklCQkKK+xs1apS89NJLfqcDAIAsEah1ru9AuvHGG01QPnnypHzyySfStWtX0x6dVkOGDJEBAwZ4lKhLlSoVoNQCAJBJ26jTSkvNFSpUMP+vVauWbNy4USZNmiTt2rWTCxcuSGJiokepWnt9Fy9ePMX9RUZGmgUAgJAM1O6l1asZP368pKXErm3MGrSzZ89u7nGtw7KUDgM7ePCgqSoHACAU+B2of/rpJ7NoRy+ttla7du2SbNmySc2aNT3arlNTTd2iRQvTQez06dOmh/eqVatk2bJlEh0dbeYQ1wuDggULmnHWffv2NUGaHt8AgFDhd6DWIVV63+lZs2ZJgQIFzDqd+KR79+7SqFEjeeaZZ1K9r+PHj5uZzLTXuAZmnfxEg3Tz5s3N8xMmTJDw8HBTotZStt61a8qUKf4mGQCA0LrN5ddff+0K0kr/P2LECHOvan8CtY6TvpKoqCiZPHmyWQAACEV+z0ymvaj/+OOPZOt1nVZfAwCAIAbq+++/31RzL1y4UA4fPmyWTz/91LQnM70nAABBrvqeNm2aDBw4UDp27Gg6lJmdRESYQP3aa68FOHkAAIQ2vwO13spSO3RpUN67d69Zd8MNN0ju3LmvRfoAAAhpab57lvbU1kXvnKVBWuf8BgAAQQ7Uf/31l9x5551SqVIladmypQnWSqu+/enxDQAArkGg1ntE64xhOkOYVoPbdMrPuLg4f3cHAAAC2UatY6h1UhK9X7Q7rQI/cOCAv7sDAACBLFHr/aLdS9K2v//+m5thAAAQYH4Hap0mdPbs2R5zeuuNNMaOHStNmjQJdPoAAAhpfld9a0DWzmSbNm0yt6EcNGiQbN++3ZSo161bd21SCQBAiPK7RH3zzTebu2U1bNhQWrdubarCdUYyvaOWjqcGAABBKlHrTGR33323mZ1s6NChAUwGAABId4lah2X9/PPP/rwEAABkZNX3I488ctXbUwIAgCB1Jvv333/l/fffl+XLl0utWrWSzfE9fvz4ACUNAACkKlBrdbd2IgsPD5dffvlFatasadZrpzJ3OlQLAABkcKCuUaOGmdO7aNGiZvaxjRs3SqFChQKYDAAAkOY26vz588v+/fvN/3/77TczwQkAAHBIibpt27bSuHFjKVGihKnerl27tmTLls3ntvv27Qt0GgEACFmpCtTvvPOOmdRkz5498tRTT0nPnj0lb9681z51AACEuFT3+taJTtTmzZvl6aefJlADAODE4VkzZsy4NikBAADpn/AEAACESKAeNWqU3HbbbaYaXYd+tWnTRnbu3Omxzfnz56V3795mOFiePHlMx7Zjx44FLc0AAIRMoF69erUJwhs2bJBvvvnG3PTjrrvuMnfksvXv31+WLFkiCxYsMNsfOXLEdGwDACAU+N1GHUhxcXEej2fOnGlK1tph7fbbb5eTJ0+aecXnzJkjTZs2dbWRV6lSxQT3evXqBSnlAACEYBu1BmZVsGBB81cDtpaymzVr5tqmcuXKUrp0aYmPj/e5j6SkJDl16pTHAgBAZuWYQK2znfXr108aNGhg5hVXCQkJkiNHDjMzmrtixYqZ51Jq946OjnYtpUqVypD0AwCQpQO1tlXrDT/mzZuXrv0MGTLElMzt5dChQwFLIwAAIdVGbevTp4988cUXsmbNGrn++utd64sXLy4XLlyQxMREj1K19vrW53yJjIw0CwAAWUFQS9SWZZkgvWjRIlm5cqWUK1fO43m933X27NllxYoVrnU6fOvgwYMSExMThBQDABBCJWqt7tYe3Z9//rkZS223O2vbcs6cOc3fHj16yIABA0wHs3z58knfvn1NkKbHNwAgFAQ1UE+dOtX8veOOOzzW6xCsbt26mf9PmDBBwsPDzUQn2qM7NjZWpkyZEpT0AgAQUoFaq76vJioqSiZPnmwWAABCjWN6fQMAgOQI1AAAOBiBGgAAByNQAwDgYARqAAAcjEANAICDEagBAHAwAjUAAA5GoAYAwMEI1AAAOBiBGgAAByNQAwDgYARqAAAcjEANAICDEagBAHAwAjUAAA5GoAYAwMEI1AAAOBiBGgAAByNQAwDgYARqAAAcjEANAICDEagBAHAwAjUAAA4W1EC9Zs0aadWqlZQsWVLCwsLks88+83jesix54YUXpESJEpIzZ05p1qyZ7N69O2jpBQAgpAL12bNn5dZbb5XJkyf7fH7s2LHyxhtvyLRp0+T777+X3LlzS2xsrJw/fz7D0woAQDBESBC1aNHCLL5oaXrixIkybNgwad26tVk3e/ZsKVasmCl5t2/fPoNTCwBAxnNsG/X+/fslISHBVHfboqOjpW7duhIfH5/i65KSkuTUqVMeCwAAmZVjA7UGaaUlaHf62H7Ol1GjRpmAbi+lSpW65mkFACDkAnVaDRkyRE6ePOlaDh06FOwkAQCQ9QJ18eLFzd9jx455rNfH9nO+REZGSr58+TwWAAAyK8cG6nLlypmAvGLFCtc6bW/W3t8xMTFBTRsAACHR6/vMmTOyZ88ejw5kW7ZskYIFC0rp0qWlX79+MmLECKlYsaIJ3M8//7wZc92mTZtgJhsAgNAI1Js2bZImTZq4Hg8YMMD87dq1q8ycOVMGDRpkxlr36tVLEhMTpWHDhhIXFydRUVFBTDUAACESqO+44w4zXjolOlvZyy+/bBYAAEKRY9uoAQAAgRoAAEejRA0AgIMRqAEAcDACNQAADkagBgDAwQjUAAA4GIEaAAAHI1ADAOBgBGoAAByMQA0AgIMRqAEAcDACNQAADkagBgDAwQjUAAA4GIEaAAAHI1ADAOBgBGoAAByMQA0AgIMRqAEAcDACNQAADkagBgDAwQjUAAA4GIEaAAAHyxSBevLkyVK2bFmJioqSunXryg8//BDsJAEAkCEcH6g//vhjGTBggAwfPlx+/PFHufXWWyU2NlaOHz8e7KQBAHDNOT5Qjx8/Xnr27Cndu3eXm266SaZNmya5cuWS999/P9hJAwDgmosQB7tw4YJs3rxZhgwZ4loXHh4uzZo1k/j4eJ+vSUpKMovt5MmT5u+pU6cClq5SeUQuFsyWrn0Uibxs0hSIfQV6f6GwLyenLRT25eS0OXVfTk6bU/cV6P3pPgIVS+z9WJZ19Y0tB/v999/1E1jr16/3WP/ss89aderU8fma4cOHm9ewkAccAxwDHAMcA+LwPDh06NBVY6GjS9RpoaVvbdO2Xb58Wf7++28pVKiQhIWFpWvf5qqsVCk5dOiQ5MuXLwCpzfrIM/KM48yZ+G0GN8+0JH369GkpWbLkVbd1dKAuXLiwZMuWTY4dO+axXh8XL17c52siIyPN4i5//vwBTZd+QQRq8uxa4zgjzzICx1nw8iw6OjrzdybLkSOH1KpVS1asWOFRQtbHMTExQU0bAAAZwdElaqXV2F27dpXatWtLnTp1ZOLEiXL27FnTCxwAgKzO8YG6Xbt28scff8gLL7wgCQkJUr16dYmLi5NixYpleFq0Sl3Hc3tXrYM84zgLLn6b5FlWPs7CtEdZhr4jAABINUe3UQMAEOoI1AAAOBiBGgAAByNQAwDgYATqdN5Sc8GCBVK5cmWz/S233CJffvmlhBp/8mz69OnSqFEjKVCggFl03vZQvG1pWm/dOm/ePDPDXps2bSTU+JtniYmJ0rt3bylRooTppVupUqWQ+336m2c6/PXGG2+UnDlzmhm4+vfvL+fPn5dQsGbNGmnVqpWZKUx/Y5999tlVX7Nq1SqpWbOmOb4qVKggM2fOvDaJC9S83FnBvHnzrBw5cljvv/++tX37dqtnz55W/vz5rWPHjvncft26dVa2bNmssWPHWjt27LCGDRtmZc+e3dq2bZsVKvzNs44dO1qTJ0+2fvrpJ+t///uf1a1bNys6Oto6fPiwFSr8zTPb/v37reuuu85q1KiR1bp1ayuU+JtnSUlJVu3ata2WLVtaa9euNXm3atUqa8uWLVao8DfPPvroIysyMtL81fxatmyZVaJECat///5WKPjyyy+toUOHWgsXLjRzcC9atOiK2+/bt8/KlSuXNWDAAHP+f/PNN008iIuLC3jaCNRu9EYfvXv3dj2+dOmSVbJkSWvUqFE+M+/hhx+27rnnHo91devWtR577DErVPibZ97+/fdfK2/evNasWbOsUJGWPNN8ql+/vvXuu+9aXbt2DblA7W+eTZ061Spfvrx14cIFK1T5m2e6bdOmTT3WaRBq0KCBFWokFYF60KBBVtWqVT3WtWvXzoqNjQ14eqj69rqlplbFpvaWmrrefXsVGxub4vZZTVryzNu5c+fk4sWLUrBgQQkFac2zl19+WYoWLSo9evSQUJOWPFu8eLGZZlirvnVypJtvvllGjhwply5dklCQljyrX7++eY1dPb5v3z7TVNCyZcsMS3dmEp+B53/Hz0yWUf7880/zI/ae8Uwf//rrrz5fozOl+dpe14eCtOSZt8GDB5s2Ie8DPqtKS56tXbtW3nvvPdmyZYuEorTkmQaZlStXSqdOnUyw2bNnjzz55JPmolBnlsrq0pJnHTt2NK9r2LChubPTv//+K48//rj897//zaBUZy4JKZz/9Q5b//zzj2nnDxRK1Aia0aNHm85RixYtMp1dkJzeBq9z586mE57eTQ6pozfv0RqId955x9zYR6ciHjp0qEybNo0svELHKK11mDJlivz444+ycOFCWbp0qbzyyivkWZBRok7HLTV1vT/bZzVpyTPb66+/bgL18uXLpVq1ahIq/M2zvXv3ym+//WZ6o7oHIRURESE7d+6UG264QbKytBxn2tM7e/bs5nW2KlWqmFKQVgvrnfmysrTk2fPPP28uCv/zn/+YxzqKRW+A1KtXL3ORo1XnuPr5X29/GcjStCLn03FLTV3vvr365ptvQuYWnGm9DenYsWPNVbreXEXvihZK/M0zHfq3bds2U+1tL/fdd580adLE/F+H0GR1aTnOGjRoYKq77YsatWvXLhPAs3qQTmueaX8R72BsX+hwS4ggn/8D3j0tkw9n0OEJM2fONN3te/XqZYYzJCQkmOc7d+5sPffccx7DsyIiIqzXX3/dDDUaPnx4SA7P8ifPRo8ebYaMfPLJJ9bRo0ddy+nTp61Q4W+eeQvFXt/+5tnBgwfNaII+ffpYO3futL744guraNGi1ogRI6xQ4W+e6flL82zu3Llm6NHXX39t3XDDDWZ0Syg4ffq0GTaqi4bG8ePHm/8fOHDAPK95pXnmPTzr2WefNed/HXbK8KwMomPhSpcubYKJDm/YsGGD67nGjRubk6S7+fPnW5UqVTLba1f9pUuXWqHGnzwrU6aM+RF4L3qSCCX+HmehHqjTkmfr1683wyU1WOlQrVdffdUMcwsl/uTZxYsXrRdffNEE56ioKKtUqVLWk08+aZ04ccIKBd9++63Pc5OdR/pX88z7NdWrVzf5q8fYjBkzrknauM0lAAAORhs1AAAORqAGAMDBCNQAADgYgRoAAAcjUAMA4GAEagAAHIxADQCAgxGoAQBwMAI1gHR58cUXpXr16uQicI0QqAEAcDACNQAADkagBrKQO+64Q/r06WOW6Ohoc19ivc+wr9sUnjp1ytw396uvvvJYv2jRIsmbN6+57aEaPHiwVKpUSXLlyiXly5c3+7t48eIV09CvXz+PdW3atJFu3bq5HiclJcnAgQPluuuuk9y5c0vdunVl1apVAcgBIOshUANZzKxZsyQiIkJ++OEHmTRpkowfP17efffdZNvpDe7vvfdemTNnjsf6jz76yARWDcxKg/bMmTNlx44dZn/Tp0+XCRMmpCuNeiERHx8v8+bNk59//lkeeughufvuu2X37t3p2i+QFUUEOwEAAqtUqVImkIaFhcmNN94o27ZtM4979uyZbNtOnTpJ586dTelZA7OWspcuXWpK1bZhw4a5/l+2bFlTEtYAO2jQoDSl7+DBgzJjxgzzt2TJkmad7jMuLs6sHzlyZJr2C2RVlKiBLKZevXomSNtiYmJMSfXVV1+VPHnyuBYNlC1btpTs2bPL4sWLzbaffvqpKWk3a9bM9fqPP/5YGjRoIMWLFzev08Ctr00rvXC4dOmSqU53T8/q1atl79696fz0QNZDiRoIEY8//ri0a9fO9VhLs1pF/uCDD5rq7/bt25u/uo2uV1o9raXul156SWJjY027t5amx40bl+L7hIeHJ2sTd2/TPnPmjGTLlk02b95s/rrTgA3AE4EayGK+//57j8cbNmyQihUrSqFChcziTQNx8+bNZfv27bJy5UoZMWKE67n169dLmTJlZOjQoa51Bw4cuOL7FylSRI4ePep6rKXnX375RZo0aWIe16hRw6w7fvy4NGrUKF2fFQgFVH0DWYxWSw8YMEB27twpc+fOlTfffFOefvrpFLe//fbbTbW2Buxy5cqZHtg2DfC6Py1Fa7X0G2+84dF+7UvTpk1NO7cuv/76qzzxxBOSmJjoel6rvPW9unTpIgsXLpT9+/ebjm+jRo0yrwHgiUANZDEaAP/55x+pU6eO9O7d2wTpXr16pbi9tmd36NBBtm7dagKou/vuu0/69+9vemnr7GNawtbhWVfy6KOPSteuXU06GjdubIZ02aVpm3Ya0+efeeYZ0+FNe5lv3LhRSpcunc5PD2Q9YZavAZYAMiUdw6wBdeLEicFOCoAAoUQNAICDEagBAHAwqr4BAHAwStQAADgYgRoAAAcjUAMA4GAEagAAHIxADQCAgxGoAQBwMAI1AAAORqAGAECc6/8BIb5OJ3FB9WYAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"mean p-value: \", pvalue_listcomp.mean())\n", "print(\"proportion < 0.05: \", (pvalue_listcomp < 0.05).mean(), \" (should be about 0.05)\")\n", "\n", "plt.figure(figsize=(5, 3.2))\n", "plt.hist(pvalue_listcomp, bins=20, edgecolor=\"white\")\n", "plt.xlabel(\"p-value\"); plt.ylabel(\"frequency\")\n", "plt.title(\"1000 null simulations: p-values are uniform\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "32691038", "metadata": {}, "source": [ "---\n", "# Cheat sheet\n", "\n", "| R | Python |\n", "|---|---|\n", "| `if (x) { } else { }` | `if x:` / `else:` (indentation, no braces) |\n", "| `else if` | `elif` |\n", "| `ifelse(cond, a, b)` | `np.where(cond, a, b)` · `a if cond else b` |\n", "| `is.na(x)` | `np.isnan(x)` · `pd.isna(x)` |\n", "| `any(z)` / `all(z)` | `any(z)` / `all(z)` |\n", "| `for (i in 1:5)` | `for i in range(1, 6):` |\n", "| `length(x)` | `len(x)` |\n", "| `x^2` | `x**2` |\n", "| `function(x) { }` | `def f(x):` |\n", "| implicit last value returned | explicit `return` required |\n", "| `sapply(x, f)` | `[f(i) for i in x]` · `list(map(f, x))` · `Series.apply(f)` |\n", "| `apply(m, 2, mean)` | `m.mean(axis=0)` |\n", "| `apply(m, 1, mean)` | `m.mean(axis=1)` |\n", "| `matrix(1:20, ncol=2)` | `np.arange(1,21).reshape(10,2, order=\"F\")` |\n", "| `rnorm(n)` | `rng.normal(size=n)` |\n", "| `sample(1:n, k)` | `rng.choice(n, size=k, replace=False)` |\n", "| `t.test(a, b)$p.value` | `stats.ttest_ind(a, b).pvalue` |\n", "| `set.seed(1)` | `np.random.default_rng(1)` |\n", "| `x[-i]` (\"drop i\") | boolean mask `x[~mask]` |\n", "| **1-based** indexing | **0-based** indexing |" ] } ], "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.13.3" } }, "nbformat": 4, "nbformat_minor": 5 }