docs_input/notebooks/03_transforms.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPU-Accelerated Numerical Computing with MatX"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial List\n",
    "1. [Introduction](01_introduction.ipynb)\n",
    "2. [Operators](02_operators.ipynb)\n",
    "3. Executors (this tutorial)\n",
    "4. [Radar Pipeline Example](04_radar_pipeline.ipynb)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Executors\n",
    "MatX executors are a generic name given to functions that execute work on the device. Operators and generators were introduced in the last tutorial as a way to generate a CUDA kernel from an expression, but they do not execute any work on the device. The `run` took an operator as input and executed it on the device. Many other types of executors exist in MatX where more complex functions can be executed alongside operators. Some executors are wrappers around existing CUDA libraries, while others are custom executors inside of MatX. This distinction is hidden from developers so that the implementation of an executor can change over time without modifying the client code. Some executors can take an operator as input, while others can only take tensors as input. These restrictions are noted in the MatX documentation, and may be relaxed or removed in future versions.\n",
    "\n",
    "Besides `run`, other executors typically allow non-element-wise kernels to execute using highly-optimized library backends. Some examples of this would be a matrix multiply (GEMM), reduction, FFT, sorting, and linear solvers. Besides the type of inputs allowed, executors may also have restrictions on the rank and/or size of a tensor. For example, performing a GEMM requires that the tensors are at least rank 2 (i.e. be a matrix), and the last dimension of the first tensor must match the second-to-last dimension of the second tensor (`MxK * KxN`). Most executors support batching, and anything above the nominal rank will result in batching dimensions. In a 1D FFT this would mean that any dimension above 1 is treated as another 1D batched FFT, and a 2D FFT would batch any dimensions above 2. \n",
    "\n",
    "Some executors use CUDA libraries to implement their functionality, and those libraries require either a handle or a plan to operated. MatX hides this complexity by creating and caching the plan on the first call, and using the same plan on future calls where possible. More advanced users may use the handle interface directly to avoid the caching. Only the caching interface will be covered in this tutorial since it's the recommended approach, but the non-cached version can be found in the documentation."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Multiply\n",
    "The `matmul` executor performs the matrix-matrix multiply of $$C = {\\alpha}A * B + {\\beta}C$$ where `A` is of dimensions `MxK`, `B` is `KxN`, and `C` is `MxN`. We first populate the `A` and `B` matrices with random values before the multiply as we did in the example above, then the GEMM is performed. Since the random number generator allocates memory sufficient to randomize the entire tensor, we create a random number generator large enough to generate values for both A or B. This allows us to create a single random number generator, but pull different random values for A and B by simply calling `run` twice. As mentioned above, any rank above 2 is consiered a batching dimension.\n",
    "\n",
    "We use rectangular matrices for `A` and `B`, while `C` will be a square matrix due to the outer dimensions of `A` and `B` matching. \n",
    "\n",
    "```c++\n",
    "(A = random<float>({8, 4}, NORMAL)).run();  \n",
    "(B = random<float>({4, 8}, NORMAL)).run();  \n",
    "\n",
    "matmul(C, A, B);\n",
    "```  \n",
    "\n",
    "Open the file [exercises/example3_gemm.cu](exercises/example3_gemm.cu) and edit the contents where you see TODO markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_gemm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "A:\n",
    "000000: -0.9247 -0.4253 -2.6438 0.1452\n",
    "000001: -0.1209 -0.5797 -0.6229 -0.3284\n",
    "000002: -1.0745 -0.3631 -1.6711 2.2655\n",
    "000003: 0.3117 -0.1842 1.2866 1.1820\n",
    "000004: -0.1271 1.2169 1.4353 1.0605\n",
    "000005: -0.4941 -1.4244 -0.7244 -1.2973\n",
    "000006: 0.0697 -0.0074 1.8969 0.6878\n",
    "000007: -0.0779 -0.8373 1.3506 -0.2879\n",
    "B:\n",
    "000000: 0.9911 1.0676 -0.6272 0.3202 -0.3110 -0.3441 -1.1709 -0.5371\n",
    "000001: 1.3390 -0.2401 1.2149 -0.2052 1.2999 0.2181 -1.2135 -1.3723\n",
    "000002: -0.4635 -0.4089 -0.0032 0.2967 -0.3587 -1.0455 -0.0450 -0.0985\n",
    "000003: 1.7608 0.9107 0.0288 -1.1128 0.0929 -0.1502 -0.9854 0.7889\n",
    "C:\n",
    "000000: -0.0050 0.3283 0.0760 -1.1547 0.6966 2.9677 1.5747 1.4554\n",
    "000001: -1.1856 -0.0342 -0.6359 0.2609 -0.5231 0.6156 1.1966 0.6628\n",
    "000002: 3.2124 1.6864 0.3035 -3.2863 0.6721 1.6973 -0.4584 3.0275\n",
    "000003: 1.5472 0.9272 -0.3894 -0.7960 -0.6881 -1.6701 -1.3640 0.8911\n",
    "000004: 2.7056 -0.0490 1.5840 -1.0446 1.2051 -1.3507 -2.4374 -0.9065\n",
    "000005: -4.3456 -1.0707 -1.4556 1.3628 -1.5586 0.8115 3.6179 1.2680\n",
    "000006: 0.3910 -0.0732 -0.0391 -0.1788 -0.6479 -2.1121 -0.8357 0.3284\n",
    "000007: -2.3314 -0.6966 -0.9810 0.8679 -1.5754 -1.5246 1.3302 0.8306\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FFT\n",
    "MatX provides an interface to do both 1D Fast Fourier Transforms (FFTs) and 2D FFTs. Any tensor above rank 1 will be batched in a 1D FFT, and any tensor above rank 2 will be batched in a 2D FFT. FFTs may either be done in-place or out-of-place by using the same or different variables for the output and inputs. Since the tensors are strongly-typed, the type of FFT (C2C, R2C, etc) is inferred by the tensor type at compile time. Similarly, the input and output size of the executor is deduced by the type of transform, and the input/output tensors must match those sizes. There's one exception to this rule, and it's when the input FFT is to be zero-padded at the end. In this case, the input tensor can be shorter than the output tensor, and the input will be zero-padded to the length of the output tensor. This is a common tactic used in signal and image processing for both speed and FFT resolution.\n",
    "\n",
    "In this example, we execute a 1D batched FFT on a 2D tensor populated with random complex floating point data. Since the FFT executor is performed in-place, the input and output types of the tensors are the same, and the type of the FFT is inferred as a complex-to-complex (`C2C`). The FFT length is specified by the inner dimension of the tensor, or 4 in this example, and the outer dimension is the number of batches, or 2. After the FFT completes, we perform on IFFT on the same tensor using the `ifft` interface. Ignoring floating point inaccuracies, the result of `ifft(fft(A))` should be the same as `A`, and this is shown by printing the tensors at each step. To perform a batched FFT on columns instead of rows, the tensor can be transposed by calling the `Permute` function used in the first tutorial. When the library detects a permuted tensor is being used, it can use technique to speed the FFT up over the naive method of converting the data in memory.\n",
    "\n",
    "```c++\n",
    "C.print();\n",
    "fft(C, C);\n",
    "C.print();\n",
    "ifft(C, C);  \n",
    "C.print();\n",
    "```\n",
    "\n",
    "Open the file [exercises/example3_1dfft.cu](exercises/example3_1dfft.cu) and edit the contents where you see TODO markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_1dfft"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "Initial C tensor:\n",
    "000000: -0.9247+0.9911j -0.4253+1.0676j -2.6438-0.6272j 0.1452+0.3202j \n",
    "000001: -0.1209-0.3110j -0.5797-0.3441j -0.6229-1.1709j -0.3284-0.5371j \n",
    "After FFT:\n",
    "000000: -3.8487+1.7517j 2.4666+2.1889j -3.2883-1.0238j 0.9718+1.0478j \n",
    "000001: -1.6518-2.3630j 0.6950+1.1112j 0.1644-0.6007j 0.3090+0.6085j \n",
    "After IFFT and normalization:\n",
    "000000: -0.9247+0.9911j -0.4253+1.0676j -2.6438-0.6272j 0.1452+0.3202j \n",
    "000001: -0.1209-0.3110j -0.5797-0.3441j -0.6229-1.1709j -0.3284-0.5371j \n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we take the same 2D tensor and perform a 2D FFT on it. Since the rank is 2, it will not be batched as in the previous example. \n",
    "\n",
    "```c++\n",
    "C.print();\n",
    "fft2(C, C);\n",
    "C.print();\n",
    "ifft2(C, C);  \n",
    "C.print();\n",
    "```\n",
    "\n",
    "As before, the results after the IFFT closely match the original `C` tensor, but with floating point error.\n",
    "\n",
    "Open the file [exercises/example3_2dfft.cu](exercises/example3_2dfft.cu) and edit the contents where you see TODO markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_2dfft"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "Intial C tensor:\n",
    "000000: -0.9247+0.9911j -0.4253+1.0676j -2.6438-0.6272j 0.1452+0.3202j \n",
    "000001: -0.1209-0.3110j -0.5797-0.3441j -0.6229-1.1709j -0.3284-0.5371j \n",
    "After FFT:\n",
    "000000: -2.0506+1.4036j -0.0405-0.0434j -2.6438-0.6272j 0.1452+0.3202j \n",
    "000001: -2.0051+2.7593j -0.4662-0.5353j -0.6229-1.1709j -0.3284-0.5371j \n",
    "After IFFT and normalization:\n",
    "000000: -1.8493+1.9823j -0.8507+2.1352j -0.6610-0.1568j 0.0363+0.0800j \n",
    "000001: -0.2417-0.6220j -1.1595-0.6882j -0.1557-0.2927j -0.0821-0.1343j \n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reductions\n",
    "A reduction operation takes multiple values and aggregates those into a smaller number of values. Most reductions take a large number of values and reduces them to a single value. Reductions are one of the most common operations perfomed on the GPU, which means they've been heavily researched and optimized for highly-parallel processors. Modern NVIDIA GPUs have special instructions for performing reductions to give even larger speedups over naive implementations. All of these details are hidden from the user and MatX automatically chooses the optimized path based on the hardware capabilities. \n",
    "\n",
    "MatX provides a set of optimized primitives to perform reductions on tensors for many common types. Reductions are supported across individual dimensions or on entire tensors, depending on the size of the output tensor. Currently supported reduction functions are `sum`, `min`, `max`,` mean`, `any`, and `all`.\n",
    "\n",
    "#### Full Reduction\n",
    "In this example we reduce an entire tensor to a single value by applying the reduction across all dimensions of the tensor. We apply the same random initialization from previous examples on a 2D tensor `A`. Note that the output tensor must be zeroed for a `sum` reduction since that value is continually added to during the reduction. Not initializing the output tensor will give undefined results since the variables are used as accumulators throughout the reduction. With the tensor initialized, we perform both a `max` and `sum` reduction across all dimensions of the tensor:\n",
    "\n",
    "```c++\n",
    "max(MD0, A);\n",
    "sum(AD0, A);\n",
    "```\n",
    "\n",
    "Open the file [exercises/example3_full_reduce.cu](exercises/example3_full_reduce.cu) and edit the contents where you see TODO markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_full_reduce"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "A:\n",
    "000000: -0.9247 -0.4253 -2.6438 0.1452 -0.1209 \n",
    "000001: -0.5797 -0.6229 -0.3284 -1.0745 -0.3631 \n",
    "000002: -1.6711 2.2655 0.3117 -0.1842 1.2866 \n",
    "000003: 1.1820 -0.1271 1.2169 1.4353 1.0605 \n",
    "Max: 2.265505\n",
    "Sum: -0.162026\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dimensional Reductions\n",
    "Reductions can also be performed across certain dimensions instead of the whole tensor. Dimensional reductions are useful in situations where each row contains data for a different user, for example, and we wish to sum up each user's data. By setting the output tensor view to a 1D tensor, independent reductions can be performed across the input tensor where each output element corresponds to a single row reduction from the input. Using the same tensor `A` from the previous example, we only change the output tensor type to be a 1D tensor instead of a scalar:\n",
    "\n",
    "```c++\n",
    "  max(MD1, A);\n",
    "  sum(AD1, A);  \n",
    "```\n",
    "\n",
    "Printing the new reduction tensors shows the reduced values across each row of the input tensor `A`.\n",
    "\n",
    "Open the file [exercises/example3_partial_reduce.cu](exercises/example3_partial_reduce.cu) and edit the contents where you see TODO markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_partial_reduce"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "A:\n",
    "000000: -0.9247 -0.4253 -2.6438 0.1452 -0.1209 \n",
    "000001: -0.5797 -0.6229 -0.3284 -1.0745 -0.3631 \n",
    "000002: -1.6711 2.2655 0.3117 -0.1842 1.2866 \n",
    "000003: 1.1820 -0.1271 1.2169 1.4353 1.0605 \n",
    "Max:\n",
    "000000: 0.1452 \n",
    "000001: -0.3284 \n",
    "000002: 2.2655 \n",
    "000003: 1.4353 \n",
    "Sum:\n",
    "000000: -3.9695 \n",
    "000001: -2.9686 \n",
    "000002: 2.0086 \n",
    "000003: 4.7676 \n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution\n",
    "MatX supports both 1D and 2D direct convolution using the `conv1d` and `conv2d` functions. FFT-based convolution can also be performed as a combination of existing primitives as a potentially faster alternative to direct convolution for large tensors. Both forms of direct convolution take in an extra mode which specifies how much of the output is saved, where `MATX_C_MODE_FULL` saves the entire filter ramp-up and down, `MATX_C_MODE_SAME` makes the input and output tensors the same size, and `MATX_C_MODE_VALID` only keeps valid samples (when the entire filter was part of the convolution). Convolution can be used to perform a rolling average of an input by making all filter values 1/N, where N is the length of the filter. In this example, we use a filter of length 3 to create a running average of the last 3 elements:\n",
    "\n",
    "```c++\n",
    "conv1d(Co, C, filt, MATX_C_MODE_FULL, 0);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_conv1d"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected output:\n",
    "```sh\n",
    "Initial C tensor:\n",
    "000000: -0.9247\n",
    "000001: -0.4253\n",
    "000002: -2.6438\n",
    "000003: 0.1452\n",
    "000004: -0.1209\n",
    "000005: -0.5797\n",
    "000006: -0.6229\n",
    "000007: -0.3284\n",
    "000008: -1.0745\n",
    "000009: -0.3631\n",
    "000010: -1.6711\n",
    "000011: 2.2655\n",
    "000012: 0.3117\n",
    "000013: -0.1842\n",
    "000014: 1.2866\n",
    "000015: 1.1820\n",
    "After conv1d:\n",
    "000000: -0.3082\n",
    "000001: -0.4500\n",
    "000002: -1.3313\n",
    "000003: -0.9747\n",
    "000004: -0.8732\n",
    "000005: -0.1851\n",
    "000006: -0.4411\n",
    "000007: -0.5103\n",
    "000008: -0.6753\n",
    "000009: -0.5887\n",
    "000010: -1.0362\n",
    "000011: 0.0771\n",
    "000012: 0.3020\n",
    "000013: 0.7977\n",
    "000014: 0.4714\n",
    "000015: 0.7615\n",
    "000016: 0.8229\n",
    "000017: 0.3940\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to a 1D convolution, a 2D convolution does the same computation over two dimensions. A tensor of at least rank 2 is needed for a 2D convolution. Below we use a filter of all ones using the `ones` operator to demonstrate the filter can also be an operator and not an existing tensor view. The result is the sum of the four values around each cell on the input:\n",
    "\n",
    "```c++\n",
    "conv2d(Co, C, filt, MATX_C_MODE_FULL, 0);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_conv2d"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last, we mentioned above that convolution can also be done in the frequency domain using FFTs. This is the preferred method for larger tensors since FFTs are much faster than direct convolutions in large sizes, and because FFT libraries are highly-optimized. FFT convolution uses more memory than direct if the inputs are not to be destroyed since it requires running an FFT on both the input signal and filter before filtering. If not done in-place, this typically requires `2N + L - 1` new elements in memory, where N is the signal length and L is the filter length. A full FFT convolution example can be found in `fft_conv.cu` in the MatX examples, but the main convolution code is shown below:\n",
    "\n",
    "```c++\n",
    "  // Perform the FFT in-place on both signal and filter\n",
    "    (sig_freq = fft(sig_time)).run();\n",
    "    (filt_freq = fft(filt_time)).run();\n",
    "\n",
    "    (sig_freq = sig_freq * filt_freq).run();\n",
    "\n",
    "    // IFFT in-place\n",
    "    (sig_freq = ifft(sig_freq)).run();\n",
    "```\n",
    "Since the expected output size of the full filtering operation is signal_len + filter_len - 1, both the filter and signal time domain inputs are shorter than the output. This would normally require a separate stage of allocating buffers of the appropriate size, zeroing them out, copying the time domain data to the buffers, and performing the  FFT. However, MatX has an API to do all of this automatically in the library using asynchronous allocations. This makes the call have a noticeable performance hit on the first call, but subsequent calls will be close to the time without allocation. To recognize that automatic padding is wanted, MatX uses the output tensor size compared to the input tensor size to determine whether to pad the input with zeros. In this case the output signal (sig_time and filt_time) are shorter than the output tensors (sig_freq and filt_freq), so it will automatically zero-pad the input.\n",
    "\n",
    "The above expression can also be combined into a single line:\n",
    "\n",
    "```c++\n",
    "(sig_freq = ifft(fft(sig_time, filtered_size) * fft(filt_time, filtered_size))).run(stream);\n",
    "```\n",
    "\n",
    "Using the convolution property $ h*x \\leftrightarrow H \\cdot X$ we simply multiply the signals element-wise after the FFT, then do an IFFT to go back to the time domain.\n",
    "\n",
    "Next, we do the same operation in the time domain using the `conv1d` function:\n",
    "\n",
    "```c++\n",
    "conv1d(time_out, sig_time, filt_time, matxConvCorrMode_t::MATX_C_MODE_FULL, 0);\n",
    "```\n",
    "\n",
    "To match the FFT results we do a full convolution to get all the samples from the filter ramp up and ramp down. However, if we wanted either valid or same mode we could slice the FFT convolution output at the appropriate places to give the same answer. Edit the file [exercises/example3_fft_conv.cu](exercises/example3_fft_conv.cu) and add the missing code where you see TODOs. After running the verification code at the bottom will check for accuracy.\n",
    "\n",
    "Expected output:\n",
    "```sh\n",
    "Verification successful\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./exercises/compile_and_run.sh example3_fft_conv"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concludes the third tutorial on MatX. In this tutorial you learned what executors are, and how they can be applied on tensor views. In the next example you will walk through an entire radar signal processing pipeline using all the primites learned up to this point. \n",
    "\n",
    "[Start Next Tutorial](04_radar_pipeline.ipynb)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.6.9 64-bit",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": ""
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}