path: root/circuitpython/extmod/ulab/docs/ulab-programming.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-10-25T21:25:53.804315Z",
     "start_time": "2020-10-25T21:25:43.765649Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__END_OF_DEFS__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Programming ulab\n",
    "\n",
    "Earlier we have seen, how `ulab`'s functions and methods can be accessed in `micropython`. This last section of the book explains, how these functions are implemented. By the end of this chapter, not only would you be able to extend `ulab`, and write your own `numpy`-compatible functions, but through a deeper understanding of the inner workings of the functions, you would also be able to see what the trade-offs are at the `python` level.\n",
    "\n",
    "\n",
    "## Code organisation\n",
    "\n",
    "As mentioned earlier, the `python` functions are organised into sub-modules at the C level. The C sub-modules can be found in `./ulab/code/`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The `ndarray` object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General comments\n",
    "\n",
    "`ndarrays` are efficient containers of numerical data of the same type (i.e., signed/unsigned chars, signed/unsigned integers or `mp_float_t`s, which, depending on the platform, are either C `float`s, or C `double`s). Beyond storing the actual data in the void pointer `*array`, the type definition has eight additional members (on top of the `base` type). Namely, the `dtype`, which tells us, how the bytes are to be interpreted. Moreover, the `itemsize`, which stores the size of a single entry in the array, `boolean`, an unsigned integer, which determines, whether the arrays is to be treated as a set of Booleans, or as numerical data, `ndim`, the number of dimensions (`uint8_t`), `len`, the length of the array (the number of entries), the shape (`*size_t`), the strides (`*int32_t`). The length is simply the product of the numbers in `shape`.\n",
    "\n",
    "The type definition is as follows:\n",
    "\n",
    "```c\n",
    "typedef struct _ndarray_obj_t {\n",
    "    mp_obj_base_t base;\n",
    "    uint8_t dtype;\n",
    "    uint8_t itemsize;\n",
    "    uint8_t boolean;\n",
    "    uint8_t ndim;\n",
    "    size_t len;\n",
    "    size_t shape[ULAB_MAX_DIMS];\n",
    "    int32_t strides[ULAB_MAX_DIMS];\n",
    "    void *array;\n",
    "} ndarray_obj_t;\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Memory layout\n",
    "\n",
    "The values of an `ndarray` are stored in a contiguous segment in the RAM. The `ndarray` can be dense, meaning that all numbers in the linear memory segment belong to a linar combination of coordinates, and it can also be sparse, i.e., some elements of the linear storage space will be skipped, when the elements of the tensor are traversed. \n",
    "\n",
    "In the RAM, the position of the item $M(n_1, n_2, ..., n_{k-1}, n_k)$ in a dense tensor of rank $k$ is given by the linear combination \n",
    "\n",
    "\\begin{equation}\n",
    "P(n_1, n_2, ..., n_{k-1}, n_k) = n_1 s_1 + n_2 s_2 + ... + n_{k-1}s_{k-1} + n_ks_k = \\sum_{i=1}^{k}n_is_i\n",
    "\\end{equation}\n",
    "where $s_i$ are the strides of the tensor, defined as \n",
    "\n",
    "\\begin{equation}\n",
    "s_i = \\prod_{j=i+1}^k l_j\n",
    "\\end{equation}\n",
    "\n",
    "where $l_j$ is length of the tensor along the $j$th axis. When the tensor is sparse (e.g., when the tensor is sliced), the strides along a particular axis will be multiplied by a non-zero integer. If this integer is different to $\\pm 1$, the linear combination above cannot access all elements in the RAM, i.e., some numbers will be skipped. Note that $|s_1| > |s_2| > ... > |s_{k-1}| > |s_k|$, even if the tensor is sparse. The statement is trivial for dense tensors, and it follows from the definition of $s_i$. For sparse tensors, a slice cannot have a step larger than the shape along that axis. But for dense tensors, $s_i/s_{i+1} = l_i$. \n",
    "\n",
    "When creating a *view*, we simply re-calculate the `strides`, and re-set the `*array` pointer.\n",
    "\n",
    "## Iterating over elements of a tensor\n",
    "\n",
    "The `shape` and `strides` members of the array tell us how we have to move our pointer, when we want to read out the numbers. For technical reasons that will become clear later, the numbers in `shape` and in `strides` are aligned to the right, and begin on the right hand side, i.e., if the number of possible dimensions is `ULAB_MAX_DIMS`, then `shape[ULAB_MAX_DIMS-1]` is the length of the last axis, `shape[ULAB_MAX_DIMS-2]` is the length of the last but one axis, and so on. If the number of actual dimensions, `ndim < ULAB_MAX_DIMS`, the first `ULAB_MAX_DIMS - ndim` entries in `shape` and `strides` will be equal to zero, but they could, in fact, be assigned any value, because these will never be accessed in an operation.\n",
    "\n",
    "With this definition of the strides, the linear combination in $P(n_1, n_2, ..., n_{k-1}, n_k)$ is a one-to-one mapping from the space of tensor coordinates, $(n_1, n_2, ..., n_{k-1}, n_k)$, and the coordinate in the linear array, $n_1s_1 + n_2s_2 + ... + n_{k-1}s_{k-1} + n_ks_k$, i.e., no two distinct sets of coordinates will result in the same position in the linear array. \n",
    "\n",
    "Since the `strides` are given in terms of bytes, when we iterate over an array, the void data pointer is usually cast to `uint8_t`, and the values are converted using the proper data type stored in `ndarray->dtype`. However, there might be cases, when it makes perfect sense to cast `*array` to a different type, in which case the `strides` have to be re-scaled by the value of `ndarray->itemsize`.\n",
    "\n",
    "### Iterating using the unwrapped loops\n",
    "\n",
    "The following macro definition is taken from [vector.h](https://github.com/v923z/micropython-ulab/blob/master/code/numpy/vector/vector.h), and demonstrates, how we can iterate over a single array in four dimensions. \n",
    "\n",
    "```c\n",
    "#define ITERATE_VECTOR(type, array, source, sarray) do {\n",
    "    size_t i=0;\n",
    "    do {\n",
    "        size_t j = 0;\n",
    "        do {\n",
    "            size_t k = 0;\n",
    "            do {\n",
    "                size_t l = 0;\n",
    "                do {\n",
    "                    *(array)++ = f(*((type *)(sarray)));\n",
    "                    (sarray) += (source)->strides[ULAB_MAX_DIMS - 1];\n",
    "                    l++;\n",
    "                } while(l < (source)->shape[ULAB_MAX_DIMS-1]);\n",
    "                (sarray) -= (source)->strides[ULAB_MAX_DIMS - 1] * (source)->shape[ULAB_MAX_DIMS-1];\n",
    "                (sarray) += (source)->strides[ULAB_MAX_DIMS - 2];\n",
    "                k++;\n",
    "            } while(k < (source)->shape[ULAB_MAX_DIMS-2]);\n",
    "            (sarray) -= (source)->strides[ULAB_MAX_DIMS - 2] * (source)->shape[ULAB_MAX_DIMS-2];\n",
    "            (sarray) += (source)->strides[ULAB_MAX_DIMS - 3];\n",
    "            j++;\n",
    "        } while(j < (source)->shape[ULAB_MAX_DIMS-3]);\n",
    "        (sarray) -= (source)->strides[ULAB_MAX_DIMS - 3] * (source)->shape[ULAB_MAX_DIMS-3];\n",
    "        (sarray) += (source)->strides[ULAB_MAX_DIMS - 4];\n",
    "        i++;\n",
    "    } while(i < (source)->shape[ULAB_MAX_DIMS-4]);\n",
    "} while(0)\n",
    "```\n",
    "\n",
    "We start with the innermost loop, the one recursing `l`. `array` is already of type `mp_float_t`, while the source array, `sarray`, has been cast to `uint8_t` in the calling function. The numbers contained in `sarray` have to be read out in the proper type dictated by `ndarray->dtype`. This is what happens in the statement `*((type *)(sarray))`, and this number is then fed into the function `f`. Vectorised mathematical functions produce *dense* arrays, and for this reason, we can simply advance the `array` pointer. \n",
    "\n",
    "The advancing of the `sarray` pointer is a bit more involving: first, in the innermost loop, we simply move forward by the amount given by the last stride, which is `(source)->strides[ULAB_MAX_DIMS - 1]`, because the `shape` and the `strides` are aligned to the right. We move the pointer as many times as given by `(source)->shape[ULAB_MAX_DIMS-1]`, which is the length of the very last axis. Hence the the structure of the loop\n",
    "\n",
    "```c\n",
    "    size_t l = 0;\n",
    "    do {\n",
    "        ...\n",
    "        l++;\n",
    "    } while(l < (source)->shape[ULAB_MAX_DIMS-1]);\n",
    "\n",
    "```\n",
    "Once we have exhausted the last axis, we have to re-wind the pointer, and advance it by an amount given by the last but one stride. Keep in mind that in the the innermost loop we moved our pointer `(source)->shape[ULAB_MAX_DIMS-1]` times by `(source)->strides[ULAB_MAX_DIMS - 1]`, i.e., we re-wind it by moving it backwards by `(source)->strides[ULAB_MAX_DIMS - 1] * (source)->shape[ULAB_MAX_DIMS-1]`. In the next step, we move forward by `(source)->strides[ULAB_MAX_DIMS - 2]`, which is the last but one stride. \n",
    "\n",
    "\n",
    "```c\n",
    "    (sarray) -= (source)->strides[ULAB_MAX_DIMS - 1] * (source)->shape[ULAB_MAX_DIMS-1];\n",
    "    (sarray) += (source)->strides[ULAB_MAX_DIMS - 2];\n",
    "\n",
    "```\n",
    "\n",
    "This pattern must be repeated for each axis of the array, and this is how we arrive at the four nested loops listed above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Re-winding arrays by means of a function\n",
    "\n",
    "\n",
    "In addition to un-wrapping the iteration loops by means of macros, there is another way of traversing all elements of a tensor: we note that, since $|s_1| > |s_2| > ... > |s_{k-1}| > |s_k|$, $P(n1, n2, ..., n_{k-1}, n_k)$ changes most slowly in the last coordinate. Hence, if we start from the very beginning, ($n_i = 0$ for all $i$), and walk along the linear RAM segment, we increment the value of $n_k$ as long as $n_k < l_k$. Once $n_k = l_k$, we have to reset $n_k$ to 0, and increment $n_{k-1}$ by one. After each such round, $n_{k-1}$ will be incremented by one, as long as $n_{k-1} < l_{k-1}$. Once $n_{k-1} = l_{k-1}$, we reset both $n_k$, and $n_{k-1}$ to 0, and increment $n_{k-2}$ by one. \n",
    "\n",
    "Rewinding the arrays in this way is implemented in the function `ndarray_rewind_array` in [ndarray.c](https://github.com/v923z/micropython-ulab/blob/master/code/ndarray.c). \n",
    "\n",
    "```c\n",
    "void ndarray_rewind_array(uint8_t ndim, uint8_t *array, size_t *shape, int32_t *strides, size_t *coords) {\n",
    "    // resets the data pointer of a single array, whenever an axis is full\n",
    "    // since we always iterate over the very last axis, we have to keep track of\n",
    "    // the last ndim-2 axes only\n",
    "    array -= shape[ULAB_MAX_DIMS - 1] * strides[ULAB_MAX_DIMS - 1];\n",
    "    array += strides[ULAB_MAX_DIMS - 2];\n",
    "    for(uint8_t i=1; i < ndim-1; i++) {\n",
    "        coords[ULAB_MAX_DIMS - 1 - i] += 1;\n",
    "        if(coords[ULAB_MAX_DIMS - 1 - i] == shape[ULAB_MAX_DIMS - 1 - i]) { // we are at a dimension boundary\n",
    "            array -= shape[ULAB_MAX_DIMS - 1 - i] * strides[ULAB_MAX_DIMS - 1 - i];\n",
    "            array += strides[ULAB_MAX_DIMS - 2 - i];\n",
    "            coords[ULAB_MAX_DIMS - 1 - i] = 0;\n",
    "            coords[ULAB_MAX_DIMS - 2 - i] += 1;\n",
    "        } else { // coordinates can change only, if the last coordinate changes\n",
    "            return;\n",
    "        }\n",
    "    }\n",
    "}\n",
    "```\n",
    "\n",
    "and the function would be called as in the snippet below. Note that the innermost loop is factored out, so that we can save the `if(...)` statement for the last axis.\n",
    "\n",
    "```c\n",
    "    size_t *coords = ndarray_new_coords(results->ndim);\n",
    "    for(size_t i=0; i < results->len/results->shape[ULAB_MAX_DIMS -1]; i++) {\n",
    "        size_t l = 0;\n",
    "        do {\n",
    "            ...\n",
    "            l++;\n",
    "        } while(l < results->shape[ULAB_MAX_DIMS - 1]);\n",
    "        ndarray_rewind_array(results->ndim, array, results->shape, strides, coords);\n",
    "    } while(0)\n",
    "\n",
    "```\n",
    "\n",
    "The advantage of this method is that the implementation is independent of the number of dimensions: the iteration requires more or less the same flash space for 2 dimensions as for 22. However, the price we have to pay for this convenience is the extra function call."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterating over two ndarrays simultaneously: broadcasting\n",
    "\n",
    "Whenever we invoke a binary operator, call a function with two arguments of `ndarray` type, or assign something to an `ndarray`, we have to iterate over two views at the same time. The task is trivial, if the two `ndarray`s in question have the same shape (but not necessarily the same set of strides), because in this case, we can still iterate in the same loop. All that happens is that we move two data pointers in sync.\n",
    "\n",
    "The problem becomes a bit more involving, when the shapes of the two `ndarray`s are not identical. For such cases, `numpy` defines so-called broadcasting, which boils down to two rules. \n",
    "\n",
    "1. The shapes in the tensor with lower rank has to be prepended with axes of size 1 till the two ranks become equal.\n",
    "2. Along all axes the two tensors should have the same size, or one of the sizes must be 1. \n",
    "\n",
    "If, after applying the first rule the second is not satisfied, the two `ndarray`s cannot be broadcast together. \n",
    "\n",
    "Now, let us suppose that we have two compatible `ndarray`s, i.e., after applying the first rule, the second is satisfied. How do we iterate over the elements in the tensors? \n",
    "\n",
    "We should recall, what exactly we do, when iterating over a single array: normally, we move the data pointer by the last stride, except, when we arrive at a dimension boundary (when the last axis is exhausted). At that point, we move the pointer by an amount dictated by the strides. And this is the key: *dictated by the strides*. Now, if we have two arrays that are originally not compatible, we define new strides for them, and use these in the iteration. With that, we are back to the case, where we had two compatible arrays. \n",
    "\n",
    "Now, let us look at the second broadcasting rule: if the two arrays have the same size, we take both `ndarray`s' strides along that axis. If, on the other hand, one of the `ndarray`s is of length 1 along one of its axes, we set the corresponding strides to 0. This will ensure that that data pointer is not moved, when we iterate over both `ndarray`s at the same time. \n",
    "\n",
    "Thus, in order to implement broadcasting, we first have to check, whether the two above-mentioned rules can be satisfied, and if so, we have to find the two new sets strides. \n",
    "\n",
    "The `ndarray_can_broadcast` function from [ndarray.c](https://github.com/v923z/micropython-ulab/blob/master/code/ndarray.c) takes two `ndarray`s, and returns `true`, if the two arrays can be broadcast together. At the same time, it also calculates new strides for the two arrays, so that they can be iterated over at the same time. \n",
    "\n",
    "```c\n",
    "bool ndarray_can_broadcast(ndarray_obj_t *lhs, ndarray_obj_t *rhs, uint8_t *ndim, size_t *shape, int32_t *lstrides, int32_t *rstrides) {\n",
    "    // returns True or False, depending on, whether the two arrays can be broadcast together\n",
    "    // numpy's broadcasting rules are as follows:\n",
    "    //\n",
    "    // 1. the two shapes are either equal\n",
    "    // 2. one of the shapes is 1\n",
    "    memset(lstrides, 0, sizeof(size_t)*ULAB_MAX_DIMS);\n",
    "    memset(rstrides, 0, sizeof(size_t)*ULAB_MAX_DIMS);\n",
    "    lstrides[ULAB_MAX_DIMS - 1] = lhs->strides[ULAB_MAX_DIMS - 1];\n",
    "    rstrides[ULAB_MAX_DIMS - 1] = rhs->strides[ULAB_MAX_DIMS - 1];\n",
    "    for(uint8_t i=ULAB_MAX_DIMS; i > 0; i--) {\n",
    "        if((lhs->shape[i-1] == rhs->shape[i-1]) || (lhs->shape[i-1] == 0) || (lhs->shape[i-1] == 1) ||\n",
    "        (rhs->shape[i-1] == 0) || (rhs->shape[i-1] == 1)) {\n",
    "            shape[i-1] = MAX(lhs->shape[i-1], rhs->shape[i-1]);\n",
    "            if(shape[i-1] > 0) (*ndim)++;\n",
    "            if(lhs->shape[i-1] < 2) {\n",
    "                lstrides[i-1] = 0;\n",
    "            } else {\n",
    "                lstrides[i-1] = lhs->strides[i-1];\n",
    "            }\n",
    "            if(rhs->shape[i-1] < 2) {\n",
    "                rstrides[i-1] = 0;\n",
    "            } else {\n",
    "                rstrides[i-1] = rhs->strides[i-1];\n",
    "            }\n",
    "        } else {\n",
    "            return false;\n",
    "        }\n",
    "    }\n",
    "    return true;\n",
    "}\n",
    "```\n",
    "\n",
    "A good example of how the function would be called can be found in [vector.c](https://github.com/v923z/micropython-ulab/blob/master/code/numpy/vector/vector.c), in the `vector_arctan2` function:\n",
    "\n",
    "```c\n",
    "mp_obj_t vector_arctan2(mp_obj_t y, mp_obj_t x) {\n",
    "    ...\n",
    "    uint8_t ndim = 0;\n",
    "    size_t *shape = m_new(size_t, ULAB_MAX_DIMS);\n",
    "    int32_t *xstrides = m_new(int32_t, ULAB_MAX_DIMS);\n",
    "    int32_t *ystrides = m_new(int32_t, ULAB_MAX_DIMS);\n",
    "    if(!ndarray_can_broadcast(ndarray_x, ndarray_y, &ndim, shape, xstrides, ystrides)) {\n",
    "        mp_raise_ValueError(translate(\"operands could not be broadcast together\"));\n",
    "        m_del(size_t, shape, ULAB_MAX_DIMS);\n",
    "        m_del(int32_t, xstrides, ULAB_MAX_DIMS);\n",
    "        m_del(int32_t, ystrides, ULAB_MAX_DIMS);\n",
    "    }\n",
    "\n",
    "    uint8_t *xarray = (uint8_t *)ndarray_x->array;\n",
    "    uint8_t *yarray = (uint8_t *)ndarray_y->array;\n",
    "    \n",
    "    ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);\n",
    "    mp_float_t *rarray = (mp_float_t *)results->array;\n",
    "    ...\n",
    "```\n",
    "\n",
    "After the new strides have been calculated, the iteration loop is identical to what we discussed in the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contracting an `ndarray`\n",
    "\n",
    "\n",
    "There are many operations that reduce the number of dimensions of an `ndarray` by 1, i.e., that remove an axis from the tensor. The drill is the same as before, with the exception that first we have to remove the `strides` and `shape` that corresponds to the axis along which we intend to contract. The `numerical_reduce_axes` function from [numerical.c](https://github.com/v923z/micropython-ulab/blob/master/code/numerical/numerical.c) does that. \n",
    "\n",
    "\n",
    "```c\n",
    "static void numerical_reduce_axes(ndarray_obj_t *ndarray, int8_t axis, size_t *shape, int32_t *strides) {\n",
    "    // removes the values corresponding to a single axis from the shape and strides array\n",
    "    uint8_t index = ULAB_MAX_DIMS - ndarray->ndim + axis;\n",
    "    if((ndarray->ndim == 1) && (axis == 0)) {\n",
    "        index = 0;\n",
    "        shape[ULAB_MAX_DIMS - 1] = 0;\n",
    "        return;\n",
    "    }\n",
    "    for(uint8_t i = ULAB_MAX_DIMS - 1; i > 0; i--) {\n",
    "        if(i > index) {\n",
    "            shape[i] = ndarray->shape[i];\n",
    "            strides[i] = ndarray->strides[i];\n",
    "        } else {\n",
    "            shape[i] = ndarray->shape[i-1];\n",
    "            strides[i] = ndarray->strides[i-1];\n",
    "        }\n",
    "    }\n",
    "}\n",
    "```\n",
    "\n",
    "Once the reduced `strides` and `shape` are known, we place the axis in question in the innermost loop, and wrap it with the loops, whose coordinates are in the `strides`, and `shape` arrays. The `RUN_STD` macro from [numerical.h](https://github.com/v923z/micropython-ulab/blob/master/code/numpy/numerical/numerical.h) is a good example. The macro is expanded in the `numerical_sum_mean_std_ndarray` function. \n",
    "\n",
    "\n",
    "```c\n",
    "static mp_obj_t numerical_sum_mean_std_ndarray(ndarray_obj_t *ndarray, mp_obj_t axis, uint8_t optype, size_t ddof) {\n",
    "    uint8_t *array = (uint8_t *)ndarray->array;\n",
    "    size_t *shape = m_new(size_t, ULAB_MAX_DIMS);\n",
    "    memset(shape, 0, sizeof(size_t)*ULAB_MAX_DIMS);\n",
    "    int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);\n",
    "    memset(strides, 0, sizeof(uint32_t)*ULAB_MAX_DIMS);\n",
    "\n",
    "    int8_t ax = mp_obj_get_int(axis);\n",
    "    if(ax < 0) ax += ndarray->ndim;\n",
    "    if((ax < 0) || (ax > ndarray->ndim - 1)) {\n",
    "        mp_raise_ValueError(translate(\"index out of range\"));\n",
    "    }\n",
    "    numerical_reduce_axes(ndarray, ax, shape, strides);\n",
    "    uint8_t index = ULAB_MAX_DIMS - ndarray->ndim + ax;\n",
    "    ndarray_obj_t *results = NULL;\n",
    "    uint8_t *rarray = NULL;\n",
    "    ...\n",
    "\n",
    "```\n",
    "Here is the macro for the three-dimensional case: \n",
    "\n",
    "```c\n",
    "#define RUN_STD(ndarray, type, array, results, r, shape, strides, index, div) do {\n",
    "    size_t k = 0;\n",
    "    do {\n",
    "        size_t l = 0;\n",
    "        do {\n",
    "            RUN_STD1((ndarray), type, (array), (results), (r), (index), (div));\n",
    "            (array) -= (ndarray)->strides[(index)] * (ndarray)->shape[(index)];\n",
    "            (array) += (strides)[ULAB_MAX_DIMS - 1];\n",
    "            l++;\n",
    "        } while(l < (shape)[ULAB_MAX_DIMS - 1]);\n",
    "        (array) -= (strides)[ULAB_MAX_DIMS - 2] * (shape)[ULAB_MAX_DIMS-2];\n",
    "        (array) += (strides)[ULAB_MAX_DIMS - 3];\n",
    "        k++;\n",
    "    } while(k < (shape)[ULAB_MAX_DIMS - 2]);\n",
    "} while(0)\n",
    "```\n",
    "In `RUN_STD`, we simply move our pointers; the calculation itself happens in the `RUN_STD1` macro below. (Note that this is the implementation of the numerically stable Welford algorithm.)\n",
    "\n",
    "```c\n",
    "#define RUN_STD1(ndarray, type, array, results, r, index, div)\n",
    "({\n",
    "    mp_float_t M, m, S = 0.0, s = 0.0;\n",
    "    M = m = *(mp_float_t *)((type *)(array));\n",
    "    for(size_t i=1; i < (ndarray)->shape[(index)]; i++) {\n",
    "        (array) += (ndarray)->strides[(index)];\n",
    "        mp_float_t value = *(mp_float_t *)((type *)(array));\n",
    "        m = M + (value - M) / (mp_float_t)i;\n",
    "        s = S + (value - M) * (value - m);\n",
    "        M = m;\n",
    "        S = s;\n",
    "    }\n",
    "    (array) += (ndarray)->strides[(index)];\n",
    "    *(r)++ = MICROPY_FLOAT_C_FUN(sqrt)((ndarray)->shape[(index)] * s / (div));\n",
    "})\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Upcasting\n",
    "\n",
    "When in an operation the `dtype`s of two arrays are different, the result's `dtype` will be decided by the following upcasting rules:\n",
    "\n",
    "1. Operations with two `ndarray`s of the same `dtype` preserve their `dtype`, even when the results overflow.\n",
    "\n",
    "2. if either of the operands is a float, the result automatically becomes a float\n",
    "\n",
    "3. otherwise\n",
    "\n",
    "    - `uint8` + `int8` => `int16`, \n",
    "    - `uint8` + `int16` => `int16`\n",
    "    - `uint8` + `uint16` => `uint16`\n",
    "    \n",
    "    - `int8` + `int16` => `int16`\n",
    "    - `int8` + `uint16` => `uint16` (in numpy, the result is a `int32`)\n",
    "\n",
    "    - `uint16` + `int16` => `float` (in numpy, the result is a `int32`)\n",
    "    \n",
    "4. When one operand of a binary operation is a generic scalar `micropython` variable, i.e., `mp_obj_int`, or `mp_obj_float`, it will be converted to a linear array of length 1, and with the smallest `dtype` that can accommodate the variable in question. After that the broadcasting rules apply, as described in the section [Iterating over two ndarrays simultaneously: broadcasting](#Iterating_over_two_ndarrays_simultaneously:_broadcasting)\n",
    "\n",
    "Upcasting is resolved in place, wherever it is required. Notable examples can be found in [ndarray_operators.c](https://github.com/v923z/micropython-ulab/blob/master/code/ndarray_operators.c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slicing and indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An `ndarray` can be indexed with three types of objects: integer scalars, slices, and another `ndarray`, whose elements are either integer scalars, or Booleans. Since slice and integer indices can be thought of as modifications of the `strides`, these indices return a view of the `ndarray`. This statement does not hold for `ndarray` indices, and therefore, the return a copy of the array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extending ulab\n",
    "\n",
    "The `user` module is disabled by default, as can be seen from the last couple of lines of [ulab.h](https://github.com/v923z/micropython-ulab/blob/master/code/ulab.h)\n",
    "\n",
    "```c\n",
    "// user-defined module\n",
    "#ifndef ULAB_USER_MODULE\n",
    "#define ULAB_USER_MODULE                (0)\n",
    "#endif\n",
    "```\n",
    "\n",
    "The module contains a very simple function, `user_dummy`, and this function is bound to the module itself. In other words, even if the module is enabled, one has to `import`:\n",
    "\n",
    "```python\n",
    "\n",
    "import ulab\n",
    "from ulab import user\n",
    "\n",
    "user.dummy_function(2.5)\n",
    "```\n",
    "which should just return 5.0. Even if `numpy`-compatibility is required (i.e., if most functions are bound at the top level to `ulab` directly), having to `import` the module has a great advantage. Namely, only the [user.h](https://github.com/v923z/micropython-ulab/blob/master/code/user/user.h) and [user.c](https://github.com/v923z/micropython-ulab/blob/master/code/user/user.c) files have to be modified, thus it should be relatively straightforward to update your local copy from [github](https://github.com/v923z/micropython-ulab/blob/master/). \n",
    "\n",
    "Now, let us see, how we can add a more meaningful function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a new ndarray\n",
    "\n",
    "In the [General comments](#General_comments) sections we have seen the type definition of an `ndarray`. This structure can be generated by means of a couple of functions listed in [ndarray.c](https://github.com/v923z/micropython-ulab/blob/master/code/ndarray.c). \n",
    "\n",
    "\n",
    "### ndarray_new_ndarray\n",
    "\n",
    "The `ndarray_new_ndarray` functions is called by all other array-generating functions. It takes the number of dimensions, `ndim`, a `uint8_t`, the `shape`, a pointer to `size_t`, the `strides`, a pointer to `int32_t`, and `dtype`, another `uint8_t` as its arguments, and returns a new array with all entries initialised to 0. \n",
    "\n",
    "Assuming that `ULAB_MAX_DIMS > 2`, a new dense array of dimension 3, of `shape` (3, 4, 5), of `strides` (1000, 200, 10), and `dtype` `uint16_t` can be generated by the following instructions\n",
    "\n",
    "```c\n",
    "size_t *shape = m_new(size_t, ULAB_MAX_DIMS);\n",
    "shape[ULAB_MAX_DIMS - 1] = 5;\n",
    "shape[ULAB_MAX_DIMS - 2] = 4;\n",
    "shape[ULAB_MAX_DIMS - 3] = 3;\n",
    "\n",
    "int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);\n",
    "strides[ULAB_MAX_DIMS - 1] = 10;\n",
    "strides[ULAB_MAX_DIMS - 2] = 200;\n",
    "strides[ULAB_MAX_DIMS - 3] = 1000;\n",
    "\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_ndarray(3, shape, strides, NDARRAY_UINT16);\n",
    "```\n",
    "\n",
    "### ndarray_new_dense_ndarray\n",
    "\n",
    "The functions simply calculates the `strides` from the `shape`, and calls `ndarray_new_ndarray`. Assuming that `ULAB_MAX_DIMS > 2`, a new dense array of dimension 3, of `shape` (3, 4, 5), and `dtype` `mp_float_t` can be generated by the following instructions\n",
    "\n",
    "```c\n",
    "size_t *shape = m_new(size_t, ULAB_MAX_DIMS);\n",
    "shape[ULAB_MAX_DIMS - 1] = 5;\n",
    "shape[ULAB_MAX_DIMS - 2] = 4;\n",
    "shape[ULAB_MAX_DIMS - 3] = 3;\n",
    "\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_dense_ndarray(3, shape, NDARRAY_FLOAT);\n",
    "```\n",
    "\n",
    "### ndarray_new_linear_array\n",
    "\n",
    "Since the dimensions of a linear array are known (1), the `ndarray_new_linear_array` takes the `length`, a `size_t`, and the `dtype`, an `uint8_t`. Internally, `ndarray_new_linear_array` generates the `shape` array, and calls `ndarray_new_dense_array` with `ndim = 1`.\n",
    "\n",
    "A linear array of length 100, and `dtype` `uint8` could be created by the function call\n",
    "\n",
    "```c\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_linear_array(100, NDARRAY_UINT8)\n",
    "```\n",
    "\n",
    "### ndarray_new_ndarray_from_tuple\n",
    "\n",
    "This function takes a `tuple`, which should hold the lengths of the axes (in other words, the `shape`), and the `dtype`, and calls internally `ndarray_new_dense_array`. A new `ndarray` can be generated by calling \n",
    "\n",
    "```c\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_ndarray_from_tuple(shape, NDARRAY_FLOAT);\n",
    "```\n",
    "where `shape` is a tuple.\n",
    "\n",
    "\n",
    "### ndarray_new_view\n",
    "\n",
    "This function crates a *view*, and takes the source, an `ndarray`, the number of dimensions, an `uint8_t`, the `shape`, a pointer to `size_t`, the `strides`, a pointer to `int32_t`, and the offset, an `int32_t` as arguments. The offset is the number of bytes by which the void `array` pointer is shifted. E.g., the `python` statement\n",
    "\n",
    "```python\n",
    "a = np.array([0, 1, 2, 3, 4, 5], dtype=uint8)\n",
    "b = a[1::2]\n",
    "```\n",
    "\n",
    "produces the array\n",
    "\n",
    "```python\n",
    "array([1, 3, 5], dtype=uint8)\n",
    "```\n",
    "which holds its data at position `x0 + 1`, if `a`'s pointer is at `x0`. In this particular case, the offset is 1. \n",
    "\n",
    "The array `b` from the example above could be generated as \n",
    "\n",
    "```c\n",
    "size_t *shape = m_new(size_t, ULAB_MAX_DIMS);\n",
    "shape[ULAB_MAX_DIMS - 1] = 3;\n",
    "\n",
    "int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);\n",
    "strides[ULAB_MAX_DIMS - 1] = 2;\n",
    "\n",
    "int32_t offset = 1;\n",
    "uint8_t ndim = 1;\n",
    "\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_view(ndarray_a, ndim, shape, strides, offset);\n",
    "```\n",
    "\n",
    "### ndarray_copy_array\n",
    "\n",
    "The `ndarray_copy_array` function can be used for copying the contents of an array. Note that the target array has to be created beforehand. E.g., a one-to-one copy can be gotten by \n",
    "\n",
    "```c\n",
    "ndarray_obj_t *new_ndarray = ndarray_new_ndarray(source->ndim, source->shape, source->strides, source->dtype);\n",
    "ndarray_copy_array(source, new_ndarray);\n",
    "\n",
    "```\n",
    "Note that the function cannot be used for forcing type conversion, i.e., the input and output types must be identical, because the function simply calls the `memcpy` function. On the other hand, the input and output `strides` do not necessarily have to be equal.\n",
    "\n",
    "### ndarray_copy_view\n",
    "\n",
    "The `ndarray_obj_t *new_ndarray = ...` instruction can be saved by calling the `ndarray_copy_view` function with the single `source` argument. \n",
    "\n",
    "\n",
    "## Accessing data in the ndarray\n",
    "\n",
    "Having seen, how arrays can be generated and copied, it is time to look at how the data in an `ndarray` can be accessed and modified. \n",
    "\n",
    "For starters, let us suppose that the object in question comes from the user (i.e., via the `micropython` interface), First, we have to acquire a pointer to the `ndarray` by calling \n",
    "\n",
    "```c\n",
    "ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(object_in);\n",
    "```\n",
    "\n",
    "If it is not clear, whether the object is an `ndarray` (e.g., if we want to write a function that can take `ndarray`s, and other iterables as its argument), we find this out by evaluating \n",
    "\n",
    "```c\n",
    "mp_obj_is_type(object_in, &ulab_ndarray_type)\n",
    "```\n",
    "which should return `true`. Once the pointer is at our disposal, we can get a pointer to the underlying numerical array as discussed earlier, i.e., \n",
    "\n",
    "```c\n",
    "uint8_t *array = (uint8_t *)ndarray->array;\n",
    "```\n",
    "\n",
    "If you need to find out the `dtype` of the array, you can get it by accessing the `dtype` member of the `ndarray`, i.e., \n",
    "\n",
    "```c\n",
    "ndarray->dtype\n",
    "```\n",
    "should be equal to `B`, `b`, `H`, `h`, or `f`. The size of a single item is stored in the `itemsize` member. This number should be equal to 1, if the `dtype` is `B`, or `b`, 2, if the `dtype` is `H`, or `h`, 4, if the `dtype` is `f`, and 8 for `d`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boilerplate\n",
    "\n",
    "In the next section, we will construct a function that generates the element-wise square of a dense array, otherwise, raises a `TypeError` exception. Dense arrays can easily be iterated over, since we do not have to care about the `shape` and the `strides`. If the array is sparse, the section [Iterating over elements of a tensor](#Iterating-over-elements-of-a-tensor) should contain hints as to how the iteration can be implemented.\n",
    "\n",
    "The function is listed under [user.c](https://github.com/v923z/micropython-ulab/tree/master/code/user/). The `user` module is bound to `ulab` in [ulab.c](https://github.com/v923z/micropython-ulab/tree/master/code/ulab.c) in the lines \n",
    "\n",
    "```c\n",
    "    #if ULAB_USER_MODULE\n",
    "        { MP_ROM_QSTR(MP_QSTR_user), MP_ROM_PTR(&ulab_user_module) },\n",
    "    #endif\n",
    "```\n",
    "which assumes that at the very end of [ulab.h](https://github.com/v923z/micropython-ulab/tree/master/code/ulab.h) the \n",
    "\n",
    "```c\n",
    "// user-defined module\n",
    "#ifndef ULAB_USER_MODULE\n",
    "#define ULAB_USER_MODULE                (1)\n",
    "#endif\n",
    "```\n",
    "constant has been set to 1. After compilation, you can call a particular `user` function in `python` by importing the module first, i.e., \n",
    "\n",
    "```python\n",
    "from ulab import numpy as np\n",
    "from ulab import user\n",
    "\n",
    "user.some_function(...)\n",
    "```\n",
    "\n",
    "This separation of user-defined functions from the rest of the code ensures that the integrity of the main module and all its functions are always preserved. Even in case of a catastrophic failure, you can exclude the `user` module, and start over.\n",
    "\n",
    "And now the function:\n",
    "\n",
    "\n",
    "```c\n",
    "static mp_obj_t user_square(mp_obj_t arg) {\n",
    "    // the function takes a single dense ndarray, and calculates the \n",
    "    // element-wise square of its entries\n",
    "    \n",
    "    // raise a TypeError exception, if the input is not an ndarray\n",
    "    if(!mp_obj_is_type(arg, &ulab_ndarray_type)) {\n",
    "        mp_raise_TypeError(translate(\"input must be an ndarray\"));\n",
    "    }\n",
    "    ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(arg);\n",
    "    \n",
    "    // make sure that the input is a dense array\n",
    "    if(!ndarray_is_dense(ndarray)) {\n",
    "        mp_raise_TypeError(translate(\"input must be a dense ndarray\"));\n",
    "    }\n",
    "    \n",
    "    // if the input is a dense array, create `results` with the same number of \n",
    "    // dimensions, shape, and dtype\n",
    "    ndarray_obj_t *results = ndarray_new_dense_ndarray(ndarray->ndim, ndarray->shape, ndarray->dtype);\n",
    "    \n",
    "    // since in a dense array the iteration over the elements is trivial, we \n",
    "    // can cast the data arrays ndarray->array and results->array to the actual type\n",
    "    if(ndarray->dtype == NDARRAY_UINT8) {\n",
    "        uint8_t *array = (uint8_t *)ndarray->array;\n",
    "        uint8_t *rarray = (uint8_t *)results->array;\n",
    "        for(size_t i=0; i < ndarray->len; i++, array++) {\n",
    "            *rarray++ = (*array) * (*array);\n",
    "        }\n",
    "    } else if(ndarray->dtype == NDARRAY_INT8) {\n",
    "        int8_t *array = (int8_t *)ndarray->array;\n",
    "        int8_t *rarray = (int8_t *)results->array;\n",
    "        for(size_t i=0; i < ndarray->len; i++, array++) {\n",
    "            *rarray++ = (*array) * (*array);\n",
    "        }\n",
    "    } else if(ndarray->dtype == NDARRAY_UINT16) {\n",
    "        uint16_t *array = (uint16_t *)ndarray->array;\n",
    "        uint16_t *rarray = (uint16_t *)results->array;\n",
    "        for(size_t i=0; i < ndarray->len; i++, array++) {\n",
    "            *rarray++ = (*array) * (*array);\n",
    "        }\n",
    "    } else if(ndarray->dtype == NDARRAY_INT16) {\n",
    "        int16_t *array = (int16_t *)ndarray->array;\n",
    "        int16_t *rarray = (int16_t *)results->array;\n",
    "        for(size_t i=0; i < ndarray->len; i++, array++) {\n",
    "            *rarray++ = (*array) * (*array);\n",
    "        }\n",
    "    } else { // if we end up here, the dtype is NDARRAY_FLOAT\n",
    "        mp_float_t *array = (mp_float_t *)ndarray->array;\n",
    "        mp_float_t *rarray = (mp_float_t *)results->array;\n",
    "        for(size_t i=0; i < ndarray->len; i++, array++) {\n",
    "            *rarray++ = (*array) * (*array);\n",
    "        }        \n",
    "    }\n",
    "    // at the end, return a micropython object\n",
    "    return MP_OBJ_FROM_PTR(results);\n",
    "}\n",
    "\n",
    "```\n",
    "\n",
    "To summarise, the steps for *implementing* a function are\n",
    "\n",
    "1. If necessary, inspect the type of the input object, which is always a `mp_obj_t` object\n",
    "2. If the input is an `ndarray_obj_t`, acquire a pointer to it by calling `ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(arg);`\n",
    "3. Create a new array, or modify the existing one; get a pointer to the data by calling `uint8_t *array = (uint8_t *)ndarray->array;`, or something equivalent\n",
    "4. Once the new data have been calculated, return a `micropython` object by calling `MP_OBJ_FROM_PTR(...)`.\n",
    "\n",
    "The listing above contains the implementation of the function, but as such, it cannot be called from `python`: \n",
    "it still has to be bound to the name space. This we do by first defining a function object in \n",
    "\n",
    "```c\n",
    "MP_DEFINE_CONST_FUN_OBJ_1(user_square_obj, user_square);\n",
    "\n",
    "```\n",
    "\n",
    "`micropython` defines a number of `MP_DEFINE_CONST_FUN_OBJ_N` macros in [obj.h](https://github.com/micropython/micropython/blob/master/py/obj.h). `N` is always the number of arguments the function takes. We had a function definition `static mp_obj_t user_square(mp_obj_t arg)`, i.e., we dealt with a single argument. \n",
    "\n",
    "Finally, we have to bind this function object in the globals table of the `user` module: \n",
    "\n",
    "```c\n",
    "STATIC const mp_rom_map_elem_t ulab_user_globals_table[] = {\n",
    "    { MP_OBJ_NEW_QSTR(MP_QSTR___name__), MP_OBJ_NEW_QSTR(MP_QSTR_user) },\n",
    "    { MP_OBJ_NEW_QSTR(MP_QSTR_square), (mp_obj_t)&user_square_obj },\n",
    "};\n",
    "```\n",
    "\n",
    "Thus, the three steps required for the definition of a user-defined function are \n",
    "\n",
    "1. The low-level implementation of the function itself\n",
    "2. The definition of a function object by calling MP_DEFINE_CONST_FUN_OBJ_N()\n",
    "3. Binding this function object to the namespace in the `ulab_user_globals_table[]`"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}