the-need-for-speed.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# The Need for Speed\n",
        "\n",
        "Joshua P. French  \n",
        "August 12, 2024\n",
        "\n",
        "To open this information in an interactive Colab notebook, click the\n",
        "Open in Colab graphic below.\n",
        "\n",
        "<a href=\"https://colab.research.google.com/github/jfrench/DataWrangleViz/blob/master/the-need-for-speed.ipynb\">\n",
        "<img src=\"https://colab.research.google.com/assets/colab-badge.svg\">\n",
        "</a>\n",
        "\n",
        "------------------------------------------------------------------------\n",
        "\n",
        "<img src=\"figures/lightning_mcqueen.jpg\">\n",
        "\n",
        "## Tips for optimizing matrix algebra operations ([YouTube](https://youtu.be/cBovzIQvw78), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=5e2bb274-0855-40ef-a60e-b09e015fec12))\n",
        "\n",
        "It is straightforward to perform matrix algebra in R. However, it is\n",
        "more difficult to write *efficient* (faster) matrix code. We provide\n",
        "some basic tips below.\n",
        "\n",
        "### `crossprod` and `tcrossprod`\n",
        "\n",
        "The `crossprod` and `tcrossprod` functions are faster approaches to\n",
        "performing matrix multiplication when one of the matrices must first be\n",
        "transposed.\n",
        "\n",
        "If we have matrices `A` and `B` of compatible dimensions, then:\n",
        "\n",
        "-   `t(A) %*% B` is equivalent to `crossprod(A, B)`.\n",
        "-   `A %*% t(B)` is equivalent to `tcrossprod(A, B)`.\n",
        "\n",
        "`crossprod` and `tcrossprod` are faster than their equivalent approaches\n",
        "because the equivalent approaches compute the transpose of matrix in\n",
        "memory before performing the multiplication. The `crossprod` and\n",
        "`tcrossprod` functions modify their multiplication algorithms so that\n",
        "the transpose matrix doesn’t have to be created prior to multiplication.\n",
        "\n",
        "To illustrate this difference, we create `A`, a $50\\times 2$ matrix of\n",
        "values and `B`, a $50 \\times 2$ matrix of values by drawing the\n",
        "appropriate number of values from a standard normal distribution. We use\n",
        "the `set.seed` function to make the example reproducible."
      ],
      "id": "8d21301a-f24b-4189-8da1-ac87ca6eb5c1"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "set.seed(314)\n",
        "A <- matrix(rnorm(100), ncol = 2)\n",
        "B <- matrix(rnorm(100), ncol = 2)"
      ],
      "id": "9f5a49d0-76c8-4e99-b89a-406631341e1d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we use the `all.equal` function to confirm that the standard and\n",
        "optimized approaches produce the same results. `all.equal` confirms that\n",
        "the supplied objects are equal to each other."
      ],
      "id": "851611ea-5a46-4eb6-b6d0-8b5abf2a5d9f"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "all.equal(t(A) %*% B, crossprod(A, B))"
      ],
      "id": "bf163738-0199-402e-a1f1-6845a9040d62"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "all.equal(A %*% t(B), tcrossprod(A, B))"
      ],
      "id": "33b70b4d-a679-46d9-8e8a-aa975554e4f4"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will use the `microbenchmark` function from the **microbenchmark**\n",
        "package (Mersmann 2023) to time the optimized and standard approaches to\n",
        "computing $\\mathbf{A}^T \\mathbf{B}$ and $\\mathbf{A} \\mathbf{B}^T$.\n",
        "\n",
        "We start by loading the **microbenchmark** package."
      ],
      "id": "9f8d6d0c-e8b2-4fa0-bba8-1b4d18ac8a8a"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(microbenchmark)"
      ],
      "id": "14785d1c-e882-461c-9d4c-296339a0bd42"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now compute $\\mathbf{A}^T \\mathbf{B}$ using the standard approach and\n",
        "using `crossprod` within the `microbenchmark` function. The\n",
        "`microbenchmark` function will perform each computation 100 times by\n",
        "default to determine the typical speed. We then use the `plot` function\n",
        "to plot the results."
      ],
      "id": "53f659de-64ee-4688-89c0-176dd5b8e7b5"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "t1 <- microbenchmark::microbenchmark(t(A) %*% B, crossprod(A, B))\n",
        "plot(t1)"
      ],
      "id": "e1fc6b74-1563-4682-b9c6-94c37c426cc8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "While the results shown in the above plot will depend on the speed of\n",
        "the computer used to perform the calculation, the boxplot shown for\n",
        "`crossprod(A, B)` should be lower then `t(A) %*% B`, which means the\n",
        "`crossprod` approach is faster than the standard approach.\n",
        "\n",
        "We perform a simialr analysis for computing $\\mathbf{A}\\mathbf{B}^T$."
      ],
      "id": "b4946cf4-7d01-40b4-95d5-11644c3352e7"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "t2 <- microbenchmark(A %*% t(B), tcrossprod(A, B))\n",
        "plot(t2)"
      ],
      "id": "25ec0009-a344-4c3b-922a-bc3af6fa893c"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once again, the plot of our timing results should show that\n",
        "`tcrossprod(A, B)` can be computed faster than `A %*% t(B)`.\n",
        "\n",
        "### Order of matrix multiplication\n",
        "\n",
        "Matrix multiplication is associative, meaning that if we have compatible\n",
        "matrices $\\mathbf{A}$, $\\mathbf{B}$ and $\\mathbf{C}$, then\n",
        "$$(\\mathbf{A B})\\mathbf{C} = \\mathbf{A} (\\mathbf{B C}).$$\n",
        "\n",
        "The order of matrix multiplication can greatly impact how long it takes\n",
        "to perform the operations.\n",
        "\n",
        "**Prioritize performing operations that produce smaller matrices\n",
        "first.** Smaller can mean the total number of rows/columns or long/thin\n",
        "matrices, depending on the context.\n",
        "\n",
        "We create three matrices below filled with random values: `A` is a\n",
        "matrix of size $500\\times 2$, `B` is a matrix of size $2 \\times 50$, `C`\n",
        "is a $50\\times 1$ column vector."
      ],
      "id": "900aef1e-3469-426b-bc02-3acf7a03ba91"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "set.seed(1999)\n",
        "A <- matrix(rnorm(1000), ncol = 2)\n",
        "B <- matrix(rnorm(100), nrow = 2)\n",
        "C <- rnorm(50)"
      ],
      "id": "8bb14c33-e611-497f-b952-61bc1b783440"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Computing $\\mathbf{A}(\\mathbf{BC})$ is much faster than computing\n",
        "$(\\mathbf{AB})\\mathbf{C}$ because $\\mathbf{BC}$ first produces a\n",
        "$2\\times 1$ matrix while $\\mathbf{AB}$ produces a much larger\n",
        "$500 \\times 50$ matrix.\n",
        "\n",
        "We confirm this by timing these operations using the `microbenchmark`\n",
        "function and plotting the results."
      ],
      "id": "e3835fb5-9739-49dd-aef0-7acb06328d4f"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "t3 <- microbenchmark((A %*% B) %*% C,\n",
        "                     A %*% (B %*% C))\n",
        "plot(t3, las = 2, cex.axis = 0.5)"
      ],
      "id": "03852235-9bc8-4bdb-ba3b-7391d6fafaff"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Loops versus vectorization ([YouTube](https://youtu.be/0AlZW5ZGvvg), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=4bf55aac-da5d-4824-92f2-b0af01266aec))\n",
        "\n",
        "Loops in R are notoriously slow.\n",
        "\n",
        "When possible, we should prefer *vectorized operations* over loops.\n",
        "\n",
        "A **vectorized** operation uses a fast, pre-compiled function to apply\n",
        "an operation to all elements of a vector, matrix, list, etc., without\n",
        "using a loop.\n",
        "\n",
        "Here is a trivial example demonstrating the speed and power of\n",
        "vectorization over loops.\n",
        "\n",
        "Let’s say we want to compute the cumulative sum of a set of numeric\n",
        "values. The cumulative sum of a vector of numeric values is the sum of\n",
        "each value and all the previous values in the vector. The cumulative sum\n",
        "of the vector $[1, 2, 4]$ is $[1, 3, 7]$.\n",
        "\n",
        "We create a function, `cumsum_loop`, to compute the cumulative sum of a\n",
        "numeric vector. We describe our approach:\n",
        "\n",
        "-   The function takes, `x`, which we assume is a numeric vector.\n",
        "-   We create a vector, `results`, to store our cumulative sums.\n",
        "-   We create a `for` loop to iterate `i` over every element of `x`. For\n",
        "    each value of `i`, we:\n",
        "    -   Check if we are considering the first element of `x`. If we are,\n",
        "        then `results[i]` is simply the value in the first position of\n",
        "        `x`.\n",
        "    -   If we are not considering the first element of `x`, then the\n",
        "        cumulative sum of the first `i` elements of `x` is the value in\n",
        "        the $i$th position of `x`, `x[i]`, plus the previous value of\n",
        "        the cumulative sum, `results[i - 1]`.\n",
        "-   After we iterating over all positions of `x`, we return the\n",
        "    `results` vector."
      ],
      "id": "4acb9054-134d-4297-9540-c120f4739b44"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum_loop <- function(x) {\n",
        "  # create vector of results\n",
        "  results <- numeric(length(x))\n",
        "  for (i in 1:length(x)) {\n",
        "    if (i == 1) {\n",
        "      results[i] <- x[i]\n",
        "    } else {\n",
        "      results[i] <- results[i - 1] + x[i]\n",
        "    }\n",
        "  }\n",
        "  return(results)\n",
        "}"
      ],
      "id": "5b8695d6-acaf-4c30-9611-f4d2fe1e729e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s try this function on the vector $[1, 2, 4]$."
      ],
      "id": "89fd5a7b-e0e6-40f7-a94f-45d547d26f87"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum_loop(c(1, 2, 4))"
      ],
      "id": "a0cb4f11-c24e-4efe-9962-eb94aea0fa15"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What is wrong with this function? Nothing, really. Except that it’s slow\n",
        "compared to the equivalent vectorized version of this calculation.\n",
        "\n",
        "The `cumsum` function is a vectorized version of a cumulative sum that\n",
        "relies on fast, pre-compiled code. We can see that it gives the same\n",
        "results as the function we created."
      ],
      "id": "99d95f0b-f12d-4f74-9309-4e23a2662212"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "cumsum(c(1, 2, 4))"
      ],
      "id": "59a348c0-7ab1-4311-a774-5247b4b50c0e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the `microbenchmark` function in the **microbenchmark** package\n",
        "(Mersmann 2023) to compare the timings of `cumsum_loop` and `cumsum`.\n",
        "The `microbenchmark` function will run the provided functions 100 times,\n",
        "by default.\n",
        "\n",
        "-   We sample 1000 values from a standard normal distribution.,\n",
        "-   We use `microbenchmark` to time how long it takes to compute the\n",
        "    cumulative sum of the sampled values using both cumulative sum\n",
        "    functions. We repeat this timing 100 times for each function.\n",
        "-   We plot the results for comparison."
      ],
      "id": "add85eb8-4820-4d47-8e22-b1de496fcd86"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# attach microbenchmark\n",
        "library(microbenchmark) \n",
        "# draw 1000 values from a N(0, 1) distribution\n",
        "x <- rnorm(1000)\n",
        "# compute the cumulative sum 100 times using both approaches\n",
        "cumsum_timings <- microbenchmark(cumsum_loop(x), cumsum(x))\n",
        "# plot timing results\n",
        "plot(cumsum_timings)"
      ],
      "id": "9a971596-1a08-4ffa-b63c-2b639cbabc88"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "While the exact timing results will depend on the computer the timings\n",
        "are performed on, the plot above will clearly show that the vectorized\n",
        "`cumsum` function is substantially faster than the `cumsum_loop`\n",
        "function.\n",
        "\n",
        "These results may seem trivial, but for large data sets where these\n",
        "operations are performed thousands or millions of times, vectorized\n",
        "operations can be many orders of magnitude faster than loops.\n",
        "\n",
        "# The vectorized \\*apply functions ([YouTube](https://youtu.be/nl-ktoJriH8), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=92d3b13f-35c1-4829-bd79-b0ad010b20ed))\n",
        "\n",
        "The `*apply` functions are a set of vectorized functions that can be\n",
        "applied different object types. The “`*`” in `*apply` is a “wildcard”\n",
        "character, meaning that it can be any set of characters.\n",
        "\n",
        "`*apply` functions will typically perform operations much faster than\n",
        "equivalent operations implemented using loops. `*apply` functions also\n",
        "make code simpler and more elegant.\n",
        "\n",
        "<a href=\"#tbl-apply\" class=\"quarto-xref\">Table 1</a> summarizes the\n",
        "various `*apply` functions available by default in R.\n",
        "\n",
        "| name | description |\n",
        "|:-----|:-----------------------------------------------------------------|\n",
        "| `apply` | Applies a function to the margins of a data frame, matrix, or array. |\n",
        "| `lapply` | Applies a function to each element of a list. |\n",
        "| `sapply` | A version of `lapply` that attempts to simplify the results to a vector, matrix, or array. |\n",
        "| `vapply` | A safer version of `sapply` that returns a pre-specified type of return value. |\n",
        "| `tapply` | Applies a function to a vector or data frame based on the levels of a categorical variable |\n",
        "| `mapply` | A multivariate version of `sapply`. |\n",
        "\n",
        "Table 1: A summary of the various `*apply` functions available in R.\n",
        "\n",
        "## The `apply` function ([YouTube](https://youtu.be/iWq_EsrXhjQ), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=98faa071-9645-4c1c-afe8-b0ad010c0e64))\n",
        "\n",
        "The `apply` function applies a function to a margin (dimension) of a\n",
        "data frame, matrix, or array.\n",
        "\n",
        "The primary arguments of the `apply` function are:\n",
        "\n",
        "-   `X`: A data frame, `matrix`, or `array`.\n",
        "-   `MARGIN`: a vector indicating the dimension to apply the function\n",
        "    to. `1` indicates rows, `2` indicates columns, `3` indicates the\n",
        "    third dimension of an array, `1:2` indicates row and columns, etc.\n",
        "-   `FUN`: the function to be applied to the margins of `X`.\n",
        "-   `...`: additional arguments that need to be passed to `FUN`.\n",
        "-   `simplify`: a logical value indicating whether results should be\n",
        "    simplified, if possible. The default is `TRUE`.\n",
        "\n",
        "We provide several examples using the `apply` function below.\n",
        "\n",
        "Consider the matrix `mat` below."
      ],
      "id": "209b9b33-59f1-4b4c-a4f8-2a02e40a8451"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "mat <- cbind(1:5, 6:10)\n",
        "mat"
      ],
      "id": "c414444b-150b-4f11-aab5-756cf6f76f81"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compute the mean value for each row of `mat` by running the\n",
        "command below."
      ],
      "id": "101c421e-e096-4e43-9406-aa2f9ed78cba"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(X = mat, MARGIN = 1, FUN = mean)"
      ],
      "id": "63681fd7-f90d-41d3-9652-d91380b3594f"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can compute the mean value for each column of `mat` by running the\n",
        "command below."
      ],
      "id": "492c923d-665e-408b-aa26-d6614d28265b"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(mat, 2, mean)"
      ],
      "id": "fd7261d5-ee02-40c6-b686-5e3dca57b1a3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that the `rowMeans` and `colMeans` functions would have produced\n",
        "the same results, with simpler code. However, the examples above are\n",
        "simple examples to illustrate usage of the `apply` function.\n",
        "\n",
        "What if we want to apply a function that has arguments we need to\n",
        "specify? E.g., let’s say we want to compute the 0.8 quantile of a set of\n",
        "values in a numeric vector `x`. We could run the command\n",
        "`quantile(x, prob = 0.8)`. How can we use a function like this with the\n",
        "`apply` function?\n",
        "\n",
        "The code below uses `apply` to compute the 0.8 quantile for each column\n",
        "of `mat`."
      ],
      "id": "c0d3fa8e-1095-4130-8c17-fd9d8d53fa2d"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(mat, 2, quantile, prob = 0.8)"
      ],
      "id": "feb4f1c8-05df-4644-a607-7c67b2b0f0ee"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, we consider a basic application of the `apply` function a\n",
        "3-dimensional array.\n",
        "\n",
        "We create a 3-dimensional array with dimension $3\\times 2 \\times 4$\n",
        "below."
      ],
      "id": "e865b44b-cae9-45a3-804b-cbd5974e3069"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "A <- array(1:24, dim = c(3, 2, 4))\n",
        "print(A)"
      ],
      "id": "d0177a78-343d-485f-83ac-d4e83e48b4b0"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s compute means for each row of `A`."
      ],
      "id": "f8729ddc-06b8-4f6d-bf35-840edbd5363c"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 1, mean)"
      ],
      "id": "229c0c04-a03d-4b95-a03a-706a731e5a21"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`11.5` comes from averaging all the values in the 1st row of `A` across\n",
        "all columns and the third dimension:\n",
        "`(1 + 4 + 7 + 10 + 13 + 16 + 19 + 22)/8 = 11.5`. More simply, it is\n",
        "equivalent to running the following command."
      ],
      "id": "6d6f4371-7325-45b0-8035-948d2ef0f604"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "mean(A[1,,])"
      ],
      "id": "f7080d65-3854-42d4-bd15-e62084fe2df8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`12.5` and `13.5` come from averaging all the values in the second and\n",
        "third rows of `A`, respectively.\n",
        "\n",
        "Let’s compute the sums across the third dimension of `A`."
      ],
      "id": "8fca6165-f800-45fa-a70c-d14949813c67"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 3, sum)"
      ],
      "id": "6d1da354-51b4-4990-923b-5ef05192cf7c"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`21` comes from adding all the values in the first element of the third\n",
        "dimension of `A`. More specifically, it is equivalent to running the\n",
        "following command."
      ],
      "id": "d42a6218-4b6b-4f38-ac52-16a581746433"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "sum(A[,,1])"
      ],
      "id": "9c5f0010-e5ce-4b02-9f52-9d71a258434b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s compute the means across the row and columns of `A`."
      ],
      "id": "dfa9b71f-1c6c-4d1b-8ebf-58ab6d3a7fe5"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "apply(A, 1:2, mean)"
      ],
      "id": "8e589b24-707f-48d7-871b-ffce43a9c94d"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The command above returns the values obtained by averaging all of the\n",
        "values in a particular row and column of `A` across the third dimension\n",
        "of `A`.\n",
        "\n",
        "`10` comes from averaging the values in `A[1, 1, 1]`, `A[1, 1, 2]`,\n",
        "`A[1, 1, 3]`, `A[1, 1, 4]`, i.e., `(1 + 7 + 13 + 19)/4 = 10`.\n",
        "Alternatively, it is equivalent to running the following command."
      ],
      "id": "e9f1cd36-fd2f-47c6-ae86-7f6fc8a92451"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "mean(A[1, 1, ])"
      ],
      "id": "5fe028c6-db47-4f6e-b3f7-60f943fbfdde"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The `lapply`, `sapply`, and `vapply` functions ([YouTube](https://youtu.be/4RjlOShPeU0), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=53c2f52e-8419-4846-83f7-b0ad010c9c25))\n",
        "\n",
        "The `lapply` function applies a function to each element of an atomic\n",
        "vector or list. It is very convenient for performing analyses in\n",
        "simulation studies. The `lapply` function always returns a list.\n",
        "\n",
        "The arguments of the `lapply` function are:\n",
        "\n",
        "-   `X`: An atomic vector or list.\n",
        "-   `FUN`: the function to be applied to each element of `X`.\n",
        "-   `...`: additional arguments that need to be passed to `FUN`.\n",
        "\n",
        "We start with a simple example that computes the sample mean of a `list`\n",
        "containing numeric vectors.\n",
        "\n",
        "We first define the our `list`."
      ],
      "id": "b421eaaf-8779-4a01-869d-08d461ee6eb4"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "simple_list <- list(1:10, 11:20)\n",
        "simple_list"
      ],
      "id": "d0bcfba7-ccb5-4539-b99b-163329833ce7"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then use `lapply` to compute the sample mean of each element of\n",
        "`sample_list`."
      ],
      "id": "25a77551-4dbe-4284-89ff-53ebf85d3267"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "lapply(simple_list, mean)"
      ],
      "id": "f9bb71ca-04e3-452a-b404-8b2f9e788e17"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `sapply` function attempts to coerce the results of `lapply` to a\n",
        "`vector`, `matrix`, or `array` (whichever is simplest) using the\n",
        "`simplify2array` function. Technically, the `sapply` function has two\n",
        "additional arguments compared to `lapply`:\n",
        "\n",
        "-   `simplify`: a logical value indicating whether the `sapply` should\n",
        "    attempt to simplify the result using `simplify2array`. The default\n",
        "    is `TRUE`.\n",
        "-   `USE.NAMES`: a logical value indicating whether names should be\n",
        "    associated with the results. The default is `TRUE`.\n",
        "\n",
        "We repeat the previous example using `sapply`."
      ],
      "id": "db50b870-a10f-462d-b741-f251b92cb6e4"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list, mean)"
      ],
      "id": "e7f7608b-3396-4875-a558-166d8e0cb333"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Instead of returning a `list` object with elements `5.5` and `15.5`, we\n",
        "get a vector with the same elements.\n",
        "\n",
        "The `sapply` function *seems* like it is the more convenient version of\n",
        "`lapply`, but this is a trap. The `sapply` function can produce\n",
        "unexpected results because it is not required to return an object of a\n",
        "specific type. The result may change depending on the input.\n",
        "\n",
        "Let’s create a second list, `simple_list2`, with two elements: the first\n",
        "element is a vector of 3 numeric values and the second element is a\n",
        "vector of 4 numeric values."
      ],
      "id": "e2e14310-ee9d-4787-81ad-4d09ae13e8a9"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "simple_list2 <- list(1:3, 4:7)\n",
        "simple_list2"
      ],
      "id": "9e54e1d8-8fe2-4989-acfe-44a7bb803041"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create a function, `add2`, that takes a `numeric` vector `x` and adds\n",
        "2 to all the elements in `x`."
      ],
      "id": "248845bd-af8c-4218-b60f-f494ab933a44"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "add2 <- function(x) {\n",
        "  x + 2\n",
        "}"
      ],
      "id": "7ac0f1c5-2549-43f5-a2ba-97bc429f8c41"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We confirm that `add2` works as intended."
      ],
      "id": "9a6f1c09-e3f7-4fe4-bf86-6902dfc917fc"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "add2(2:3)"
      ],
      "id": "46ba36d8-ef2d-42a0-996b-2879e6eca5fc"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we use `sapply` to apply `add2` to each element of `simple_list`, the\n",
        "result is returned as a `matrix` because the results can be simplified\n",
        "by `simplify2array`."
      ],
      "id": "cd1c459c-6749-4450-a9b7-13219fa3aa44"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list, add2)"
      ],
      "id": "2a0cf5d9-00b1-44b8-b271-598d287027b3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If we apply `add2` to each element of `simple_list2`, the result is\n",
        "returned as a `list` because the results cannot be simplified by\n",
        "`simplify2array`."
      ],
      "id": "869b148f-2e38-4739-bc6a-337234d39878"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "sapply(simple_list2, add2)"
      ],
      "id": "924a85de-e37a-494a-9c3b-866f70f015c2"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We don’t always know that data type the `sapply` function will return,\n",
        "which is dangerous if we use those results in any future calculations.\n",
        "\n",
        "The `vapply` function attempts to be a safer version of `sapply` by\n",
        "requiring the user to indicate the data type expected from each\n",
        "application of `FUN`. However, `vapply` also attempts to simplify the\n",
        "output using `simplify2array`, so it can also create issues in later\n",
        "code.\n",
        "\n",
        "The `vapply` function essentially adds one additional argument to\n",
        "`sapply`.\n",
        "\n",
        "-   `FUN.VALUE`: a template for the value returned by `FUN`.\n",
        "\n",
        "If `FUN` doesn’t have the same structure as `FUN.VALUE`, then `vapply`\n",
        "returns an error.\n",
        "\n",
        "Let’s create a trivial example. We create two lists, `l1` and `l2`.\n",
        "\n",
        "-   `l1` has two elements containing `integer` vectors.\n",
        "-   `l2` has two elements containing `double` vectors."
      ],
      "id": "0f713f1b-ff67-477e-9295-0c9d8a684678"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "l1 <- list(seq.int(1, 4), seq.int(4, 7))\n",
        "l1"
      ],
      "id": "2a67c498-1ed6-4014-83de-f8a0c831da48"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "l2 <- list(seq(1.5, 4.5, by = 1), seq(4.5, 7.5, by = 1))\n",
        "l2"
      ],
      "id": "9d314e5b-e0eb-4aad-96c3-00bbe6262fc8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We apply the `max` function to `l1` using `vapply`. We specify\n",
        "`FUN.VALUE = integer(1)` to indicate that each time we apply the `max`\n",
        "function to an element of `l1`, we expect to get an `integer` vector of\n",
        "length 1."
      ],
      "id": "a7c3b1b9-2a06-46d3-88de-22aa022b1b69"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l1, max, FUN.VALUE = integer(1))"
      ],
      "id": "170a5574-13f4-4060-b122-5719435cb17b"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We get an error if we do the same thing for `l2` because the `max`\n",
        "function will return an atomic vector of type `double` of length 1. The\n",
        "error produced by `vapply` indicates that it was expecting to get an\n",
        "`integer` type by got a `double` type instead."
      ],
      "id": "7b26ce44-6e56-4a03-90c7-224120e0bf68"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l2, max, FUN.VALUE = integer(1))"
      ],
      "id": "510ca90e-a4e2-4bfa-80c4-f4d3227493df"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can correct the problem by changing `FUN.VALUE` to `double(1)`."
      ],
      "id": "c8d17d9b-8b00-4ccc-96e7-869e84fe2ae6"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "vapply(l2, max, FUN.VALUE = double(1))"
      ],
      "id": "1c4d7ddb-45cf-4769-bd79-e927b9174aef"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We should not use `sapply` or `vapply` inside any function to avoid\n",
        "problems with the output change data types. Instead, use the `lapply`\n",
        "function to keep the output consistent.\n",
        "\n",
        "### Bootstrap confidence interval example ([YouTube](https://youtu.be/IBuJSRu9Quw), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=88db6193-d823-4c75-a00b-b0ad010d3fb3))\n",
        "\n",
        "A confidence interval procedure produces plausible values for a target\n",
        "parameter with an associated confidence level.\n",
        "\n",
        "-   If we repeatedly obtain samples from a population, a confidence\n",
        "    interval procedure will produce intervals containing the target\n",
        "    parameter with probability equal to the confidence level.\n",
        "\n",
        "Confidence intervals are typically constructed on the assumption that we\n",
        "know the population distribution from which the sample came.\n",
        "\n",
        "A bootstrap confidence interval procedure can be used to construct\n",
        "confidence intervals when the population distribution is unknown.\n",
        "\n",
        "More generally, the bootstrap method allows us to approximate the\n",
        "sampling distribution of a statistic by using the observed data to\n",
        "produce simulated data sets.\n",
        "\n",
        "Suppose we have a sample of $n$ data values, $x_1, x_2, \\ldots, x_n$.\n",
        "Suppose we want to estimate the mean, $\\mu$, of the population from\n",
        "which the sample came. A natural estimator of $\\mu$ is the sample mean,\n",
        "$\\bar{x}$.\n",
        "\n",
        "The procedure for producing a bootstrap confidence interval for $\\mu$\n",
        "with confidence level $1-\\alpha$ is:\n",
        "\n",
        "1.  Sample with replacement a data set of size $n$ from the observed\n",
        "    data values. This is a bootstrap sample.\n",
        "2.  Compute $\\bar{x}$ for the random sample.\n",
        "3.  Perform steps 1 and 2 $B$ times, where $B$ is a number like 1000 or\n",
        "    10000.\n",
        "4.  Compute the $\\alpha/2$ and $1-\\alpha/2$ quantiles of the sample\n",
        "    means obtained from the bootstrap samples.\n",
        "\n",
        "Suppose we have observed the following sample:\n",
        "\n",
        "3.318, 4.940, 5.516, 2.958, 6.098, 6.858, -0.539, 2.544, 4.091, 1.884\n",
        "\n",
        "Let’s construct a 95% bootstrap confidence interval for the population\n",
        "mean, $\\mu$, using the sample data. We will build a function to\n",
        "construct our bootstrap confidence interval that is based on `lapply`.\n",
        "Our function will take three arguments:\n",
        "\n",
        "-   `x`: A vector of numeric values.\n",
        "-   `conf.level`: The confidence level desired for the interval. The\n",
        "    default value will be `0.95`.\n",
        "-   `B`: The number of bootstrap samples from which to build the\n",
        "    confidence interval. The default value will be `1000`.\n",
        "\n",
        "The function will perform the following actions:\n",
        "\n",
        "-   Determine `alpha` from `confidence.level`.\n",
        "-   The `lapply` function is applied to each element of the atomic\n",
        "    vector `1:B`.\n",
        "-   We supply an *anonymous* function to the `FUN` argument of `lapply`.\n",
        "    An anonymous function is a temporary function we create to perform\n",
        "    an action but for which we don’t want to go to the trouble of\n",
        "    storing in memory for long-term use. In our case, the anonymous\n",
        "    function takes the argument `i`, which comes from the elements of\n",
        "    `1:B`.\n",
        "-   Inside our anonymous function, we sample `n` values with replacement\n",
        "    from `x` using the `sample` function, making sure to set the\n",
        "    `replace` argument of `sample` to `TRUE`. This produces a bootstrap\n",
        "    sample.\n",
        "-   We compute the sample mean of our bootstrap sample and return that\n",
        "    value from our anonymous function.\n",
        "-   We store the $B$ sample means produced by the `lapply` function in\n",
        "    the `bstats` object.\n",
        "-   We use the `unlist` function to convert the list of sample means in\n",
        "    `bstats` to a `numeric` vector.\n",
        "-   We use the `quantile` function to return the `alpha/2` and\n",
        "    `1-alpha/2` quantiles from our vector of sample means."
      ],
      "id": "22d6c02f-57c9-4fa9-b921-70e0a8a06b4c"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "bootstrap_ci <- function(x, conf.level = 0.95, B = 1000) {\n",
        "  # determine alpha\n",
        "  alpha <- 1 - conf.level\n",
        "  # compute sample mean of B bootstrap samples\n",
        "  bstats <- lapply(1:B, FUN = function(i) {\n",
        "    # draw bootstrap sample\n",
        "    boot_sample <- sample(x, size = length(x), replace = TRUE)\n",
        "    # compute sample mean of bootstrap sample\n",
        "    mean(boot_sample)\n",
        "  })\n",
        "  # convert list to vector\n",
        "  bstats <- unlist(bstats)\n",
        "  # compute desired quantiles for bootstrap ci\n",
        "  quantile(bstats, prob = c(alpha/2, 1-alpha/2))\n",
        "}"
      ],
      "id": "d2e4e35a-38b8-4d5a-90bc-457947a31dcd"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now apply our function to the observed sample of data."
      ],
      "id": "6c527793-3e49-4704-8e4f-8f293bcad2ea"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "set.seed(43)\n",
        "x <- c(3.318, 4.940, 5.516, 2.958, 6.098, 6.858, -0.539, 2.544, 4.091, 1.884)\n",
        "bootstrap_ci(x)"
      ],
      "id": "7f8b73be-6077-4370-b218-a25e236a0811"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Our 95% bootstrap confidence interval for the population mean is \\[2.44,\n",
        "4.98\\]. The sample was a random sample from a normal distribution with\n",
        "mean 3 and standard deviation 2, so our interval (in this case) did\n",
        "contain the true value.\n",
        "\n",
        "## The `tapply` function ([YouTube](https://youtu.be/CLLK9a8d_jA), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=1058278d-dbb8-4edd-b2ef-b0ad010de1b8))\n",
        "\n",
        "The `tapply` function is used to apply a function to all the values\n",
        "associated with a certain level of a `factor` variable. It is useful for\n",
        "summarizing information for different groups of data.\n",
        "\n",
        "The `tapply` function takes the following arguments:\n",
        "\n",
        "-   `X`: an object for which a `split` method exists. Usually, this is a\n",
        "    `vector` or `data.frame` object.\n",
        "-   `INDEX`: a `factor` variable with the same length as `X`. This will\n",
        "    be used to “split” `X` into different groups.\n",
        "-   `FUN`: The function to be applied to each group of `X` after the\n",
        "    `split` function is used.\n",
        "-   `...`: additional arguments that will be passed to `FUN`.\n",
        "-   `default`: an argument we can probably ignore. Run `?tapply` in the\n",
        "    Console for more details about this argument.\n",
        "-   `simplify`: a logical value indicating whether the result should be\n",
        "    simplified to a scalar array. The default is `TRUE`. Otherwise, a\n",
        "    `list` is always returned.\n",
        "\n",
        "We will use the `tapply` function to summarize chicken growth across\n",
        "different feed types.\n",
        "\n",
        "The `chickwts` data set in the **datasets** package is a data frame with\n",
        "71 observations of 2 variables. The data provide the results of an\n",
        "experiment to measure and compare the effectiveness of various feed\n",
        "supplements on the growth rate of chickens. The following variables are\n",
        "included in the data frame:\n",
        "\n",
        "-   `weight`: the chick’s weight (g).\n",
        "-   `feed`: a `factor` indicating the feed type.\n",
        "\n",
        "We use `tapply` to compute the average chicken weight of the chicken’s\n",
        "receiving each feed type."
      ],
      "id": "33e777b9-c714-4711-979b-a6019c8cc1e5"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "tapply(X = chickwts$weight, INDEX = chickwts$feed, FUN = mean)"
      ],
      "id": "37c0f996-5257-4dc2-99d7-454a94ac2be3"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The average weight of the chickens eating the `casein` feed type was\n",
        "323.6 g. The average weight of the chickens eating the `soybean` feed\n",
        "type was 246.4 g.\n",
        "\n",
        "We can use the `tapply` function to compute a more complex summary of\n",
        "each group’s weights using the `summary` function."
      ],
      "id": "049d9890-7b6f-4d60-80fa-b2b60058fd41"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "tapply(X = chickwts$weight, INDEX = chickwts$feed, FUN = summary)"
      ],
      "id": "bcdec2e2-6844-4ddf-8c82-1c2a529863b2"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Parallel computing, progress bars, `*apply` functions ([YouTube](https://youtu.be/2LXW8d-lOEM), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=92a4020b-f79c-456d-a91c-b0ad0110622e))\n",
        "\n",
        "We can sometimes speed up the computations performed by an `*apply`\n",
        "function using parallel processing (processing multiple tasks *in\n",
        "parallel* instead of one-at-a-time).\n",
        "\n",
        "Parallel processing is particularly helpful with using `lapply`,\n",
        "`sapply`, or `vapply` to perform complex analysis on a sequence of\n",
        "simulated data sets. The analysis can be performed on each simulated\n",
        "data set separately, which means the operation can be parallelized.\n",
        "\n",
        "We provide small examples to illustrate the advantages of parallel\n",
        "computing when using `*apply` functions.\n",
        "\n",
        "We start by creating a simple function, `sleep_iterator`.\n",
        "\n",
        "-   The function takes the argument `i`, which is not used in the\n",
        "    function, so that it is compatible with the `lapply` function.\n",
        "-   Inside the function, we use `Sys.sleep(1)` to pause for one second.\n",
        "-   We have the function return `NULL`.\n",
        "\n",
        "We define `sleep_iterator` below."
      ],
      "id": "64042763-d985-4b64-b639-6ae1399f0f3c"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "sleep_iterator <- function(i) {\n",
        "  Sys.sleep(1) # sleep for 1 second\n",
        "  return(NULL) # return NULL\n",
        "}"
      ],
      "id": "6e5781a3-77bd-448e-8669-e543e1c1e986"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we use the `lapply` function to run `sleep_iterator` 6 times.\n",
        "\n",
        "-   `lapply` will run `sleep_iterator` for `i = 1`, `2`, `...`, `6`.\n",
        "-   We time the operations performed by `lapply` using `system.time`.\n",
        "\n",
        "Unsurprisingly, applying `sleep_iterator` six times results in an\n",
        "elapsed time of about 6 seconds since we pause for 1 second for each of\n",
        "the six values of `i`."
      ],
      "id": "0dc60e8f-f984-4f8c-b6c4-7f1c1810e6c4"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(lapply(1:6, sleep_iterator))"
      ],
      "id": "2d3c2577-bf4b-4695-8337-d020f67b8b46"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The **parallel** package (R Core Team 2023) provides parallelized\n",
        "equivalents of several `*apply` functions.\n",
        "\n",
        "-   `parApply` is a parallelized version of `apply`.\n",
        "-   `parLapply` is a parallelized version of `lapply`.\n",
        "-   `parSapply` is a parallelized version of `sapply`.\n",
        "\n",
        "Run `?parallel::clusterApply` to see additional details about these and\n",
        "related parallelized functions.\n",
        "\n",
        "We load the **parallel** package."
      ],
      "id": "cfff639d-e8ce-4c3c-9f91-5514540dc192"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(parallel)"
      ],
      "id": "7789c99a-dbdb-4588-8286-ea725b97e5d9"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The available options and behavior for performing parallel processing\n",
        "depends on our computer’s operating system. In general, it is easier to\n",
        "perform parallel computing on Macs and Linux computers. Windows\n",
        "computers cannot run forked processes through R.\n",
        "\n",
        "We want to use `parLapply` to run `sleep_iterator` 6 times using\n",
        "parallel processing.\n",
        "\n",
        "The primary arguments of `parLapply` are:\n",
        "\n",
        "-   `cl`: a cluster object returned by the `makeCluster` function that\n",
        "    defines how the processes will be run in parallel.\n",
        "-   `X`: the atomic vector or list to be iterated over.\n",
        "-   `fun`: the function to be applied to each element of `X`.\n",
        "-   `...`: additional arguments to be passed to `fun`.\n",
        "\n",
        "We use `makeCluster` to setup a parallel processing approach that runs 2\n",
        "processes in parallel. The code is supposed to be general enough that it\n",
        "will work on Windows, Mac, or Linux computers."
      ],
      "id": "ba6c43dd-edf0-4fc4-9dd3-5f5c28253be4"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "cl <- makeCluster(getOption(\"cl.cores\", 2))"
      ],
      "id": "451702d7-c17c-4d91-adb5-90fbcece8cf8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We once again use `system.time` to time how long it takes to run\n",
        "`sleep_iterator` 6 times."
      ],
      "id": "9fea0090-62e8-49e6-ae55-5419a57c4ee3"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(parLapply(cl, 1:6, sleep_iterator))"
      ],
      "id": "7881a9f8-5146-423b-bdb4-81fe222bcf91"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since we are running `sleep_iterator` in parallel using 2 processes, it\n",
        "only takes about 3 seconds to complete our computations instead of the\n",
        "original 6 seconds!\n",
        "\n",
        "The **future.apply** package (Bengtsson 2023b) provides parallelized\n",
        "implementations of several `*apply` functions. The advantage of the\n",
        "**future.apply** package is that code can be created that generalizes\n",
        "from being “evaluated on the local machine, in parallel a set of local\n",
        "machines, or distributed on a mix of local and remote machines”\n",
        "Bengtsson (2021).\n",
        "\n",
        "We start by loading the **future.apply** package."
      ],
      "id": "ef308de9-cfe7-49bb-8efa-a8f580e7ffee"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(future.apply)"
      ],
      "id": "53e6094f-dbc2-4be5-af1c-39e284231e7e"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To perform a parallelized operation using **future.apply**, we must use\n",
        "the `plan` function to describe how we want the parallel processing to\n",
        "be performed.\n",
        "\n",
        "The code below sets up a plan to run 2 processes in parallel using a\n",
        "`\"multisession\"` approach. Run `?future::plan` in the Console for more\n",
        "details."
      ],
      "id": "7517d7e9-2d1d-4ced-aa8f-00f2a2619135"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "plan(multisession, workers = 2)"
      ],
      "id": "603265df-2dc1-44d9-af33-38db47eea0fb"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we use the `future_lapply` function to run the `sleep_iterator`\n",
        "function 6 times. The function works identically to the `lapply`\n",
        "function, but uses the parallelization plan to speed up the\n",
        "computations. Similar to before, it only takes about 3 seconds to run\n",
        "the `sleep_iterator` function 6 times because we are running two\n",
        "processes in parallel."
      ],
      "id": "05b91df4-92ee-4f95-a9cc-ce9c35f4381b"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(future_lapply(1:6, FUN = sleep_iterator))"
      ],
      "id": "5d163aef-1696-4086-a814-783af4d11da8"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "When an `*apply` function takes a long time to return its results, it\n",
        "can be nice to have a progress bar to assess how long the computations\n",
        "will take.\n",
        "\n",
        "The **pbapply** package (Solymos and Zawadzki 2023) creates equivalents\n",
        "of many `*apply` functions that add progress bars showing the time\n",
        "remaining to complete the computation.\n",
        "\n",
        "-   Many of these functions are equipped to use parallel processing\n",
        "    through the **parallel** package. We must set the `cl` argument with\n",
        "    a `makeCluster` object to use this functionality.\n",
        "-   Run `?pbapply::pbapply` in the Console for a complete list of\n",
        "    available functions.\n",
        "\n",
        "We start by loading the **pbapply** package."
      ],
      "id": "7a5300bd-b786-404b-9e1c-66fb9fe5d13b"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "library(pbapply)"
      ],
      "id": "65318d69-d645-4398-a198-e6aea7229743"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `pblapply` function is an equivalent of `lapply` that includes a\n",
        "progress bar. We can provide a cluster object from the `makeCluster`\n",
        "function to perform the computations in parallel.\n",
        "\n",
        "We use `pblapply` to run `sleep_iterator` 6 times. We use the previously\n",
        "created cluster object, `cl`, which runs 2 processes in parallel.\n",
        "Similar to before, our code runs in about 3 seconds while also showing a\n",
        "progress bar. (We need to run this code interactively to see this\n",
        "progress bar.)"
      ],
      "id": "91209d4f-c138-4287-b50b-43986665c576"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(pblapply(1:6, FUN = sleep_iterator, cl = cl))"
      ],
      "id": "f0633763-2ad1-4ed6-981a-eade59f61165"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let’s a non-trivial example of parallel computing where we compare the\n",
        "execution time of `bootstrap_ci` to a parallelized version\n",
        "(`bootstrap_ci_cl`)."
      ],
      "id": "322b273e-24d4-4b91-8443-1da4135e1fd5"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "bootstrap_ci_cl <- function(x, conf.level = 0.95, B = 1000, cl = NULL) {\n",
        "  # determine alpha\n",
        "  alpha <- 1 - conf.level\n",
        "  # compute sample mean of B bootstrap samples\n",
        "  bstats <- pbapply::pblapply(1:B, FUN = function(i, xvec) {\n",
        "    # draw bootstrap sample\n",
        "    boot_sample <- sample(x, size = length(x), replace = TRUE)\n",
        "    # compute sample mean of bootstrap sample\n",
        "    mean(boot_sample)\n",
        "  }, xvec = x, cl = cl)\n",
        "  # convert list to vector\n",
        "  bstats <- unlist(bstats)\n",
        "  # compute desired quantiles for bootstrap ci\n",
        "  quantile(bstats, prob = c(alpha/2, 1-alpha/2))\n",
        "}"
      ],
      "id": "153fb639-7aea-4f7b-b772-bf36b4be5f55"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We run the function without parallelization."
      ],
      "id": "e877583d-f73e-40ef-b050-c8648275bb99"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(bootstrap_ci_cl(x, B = 100000, cl = NULL))"
      ],
      "id": "b17e826a-89ba-4b37-80a8-a3dd04482604"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We run the function with parallelization."
      ],
      "id": "32fdb4e1-8e61-4999-82f7-3914fa83268d"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "system.time(bootstrap_ci_cl(x, B = 100000, cl = cl))"
      ],
      "id": "161c02d1-99e9-41ff-b6b2-0eb7fc7bb365"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Bengtsson, Henrik. 2021. “A Unifying Framework for Parallel and\n",
        "Distributed Processing in r Using Futures.” *The R Journal* 13 (2):\n",
        "208–27. <https://doi.org/10.32614/RJ-2021-048>.\n",
        "\n",
        "———. 2023a. *Future: Unified Parallel and Distributed Processing in r\n",
        "for Everyone*. <https://CRAN.R-project.org/package=future>.\n",
        "\n",
        "———. 2023b. *Future.apply: Apply Function to Elements in Parallel Using\n",
        "Futures*. <https://CRAN.R-project.org/package=future.apply>.\n",
        "\n",
        "Mersmann, Olaf. 2023. *Microbenchmark: Accurate Timing Functions*.\n",
        "<https://github.com/joshuaulrich/microbenchmark/>.\n",
        "\n",
        "R Core Team. 2023. *R: A Language and Environment for Statistical\n",
        "Computing*. Vienna, Austria: R Foundation for Statistical Computing.\n",
        "<https://www.R-project.org/>.\n",
        "\n",
        "Solymos, Peter, and Zygmunt Zawadzki. 2023. *Pbapply: Adding Progress\n",
        "Bar to \\*Apply Functions*. <https://github.com/psolymos/pbapply>."
      ],
      "id": "763eefb5-1d1e-47be-81ae-86a0ffce08ae"
    }
  ],
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  }
}