.. |spacingstart| raw:: latex

   \begin{spacing}{1.5}



.. |spacingend| raw:: latex

   \end{spacing}







.. |newpage| raw:: latex

   \newpage


.. |begin_flushleft| raw:: latex

   \begin{flushleft}


.. |end_flushleft| raw:: latex

   \end{flushleft}


.. |vspace| raw:: html

   <br />







|newpage|

Numpy indexing
==========================================================


https://numpy.org/doc/stable/user/basics.indexing.html


ndarrays can be indexed using the standard Python x[obj] syntax, where x is the array and obj the selection. There are different kinds of indexing available depending on obj: basic indexing, advanced indexing and field access.

Most of the following examples show the use of indexing when referencing data in an array. The examples work just as well when assigning to an array. See Assigning values to indexed arrays for specific examples and explanations on how assignments work.

Note that in Python, x[(exp1, exp2, ..., expN)] is equivalent to x[exp1, exp2, ..., expN]; the latter is just syntactic sugar for the former.






Single element indexing
----------------------------------------------------------------

See https://numpy.org/doc/stable/user/basics.indexing.html#basic-indexing


Single element indexing works exactly like that for other standard Python sequences. It is 0-based, and accepts negative indices for indexing from the end of the array.


.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm
    >>> x = np.arange(10)
    >>> x[2]
    2
    >>> x[-2]
    8

It is not necessary to separate each dimension’s index into its own set of square brackets.

.. code-block:: pycon

    >>> x.shape = (2, 5)  # now x is 2-dimensional
    >>> x[1, 3]
    8
    >>> x[1, -1]
    9

Note that if one indexes a multidimensional array with fewer indices than dimensions, one gets a subdimensional array. For example:


.. code-block:: pycon

    >>> x[0]
    array([0, 1, 2, 3, 4])


That is, each index specified selects the array corresponding to the rest of the dimensions selected. In the above example, choosing 0 means that the remaining dimension of length 5 is being left unspecified, and that what is returned is an array of that dimensionality and size. It must be noted that the returned array is a view, i.e., it is not a copy of the original, but points to the same values in memory as does the original array. In this case, the 1-D array at the first position (0) is returned. So using a single index on the returned array, results in a single element being returned. That is:


.. code-block:: pycon

    >>> x[0][2]
    2




Slicing and striding
----------------------------------------------------------------


Basic slicing extends Python’s basic concept of slicing to N dimensions. Basic slicing occurs when obj is a slice object (constructed by start:stop:step notation inside of brackets), an integer, or a tuple of slice objects and integers. Ellipsis and newaxis objects can be interspersed with these as well.

The simplest case of indexing with N integers returns an array scalar representing the corresponding item. As in Python, all indices are zero-based: for the i-th index `n_i`, the valid range is `0 \le n_i < d_i`, where `d_i` is the i-th element of the shape of the array. Negative indices are interpreted as counting from the end of the array (i.e., if `n_i < 0`, it means `n_i + d_i`).

All arrays generated by basic slicing are always views of the original array. NumPy slicing creates a view instead of a copy as in the case of built-in Python sequences such as string, tuple and list. Care must be taken when extracting a small portion from a large array which becomes useless after the extraction, because the small portion extracted contains a reference to the large original array whose memory will not be released until all arrays derived from it are garbage-collected. In such cases an explicit copy() is recommended.

The standard rules of sequence slicing apply to basic slicing on a per-dimension basis (including using a step index). Some useful concepts to remember include:

The basic slice syntax is i:j:k where `i` is the starting index, `j` is the stopping index, and `k` is the step (`k \ne 0`). This selects the `m` elements (in the corresponding dimension) with index values `i, i + k, …, i + (m - 1) k` where `m = q + (r \ne 0)` and `q` and `r` are the quotient and remainder obtained by dividing `j - i` by `k: j - i = q k + r`, so that `i + (m - 1) k < j`. For example:



.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm
    >>> x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    >>> x[1:7:2]
    array([1, 3, 5])

Negative `i` and `j` are interpreted as `n + i` and `n + j` where `n` is the number of elements in the corresponding dimension. Negative `k` makes stepping go towards smaller indices. From the above example:

.. code-block:: pycon

    >>> x[-2:10]
    array([8, 9])
    >>> x[-3:3:-1]
    array([7, 6, 5, 4])

Assume `n` is the number of elements in the dimension being sliced. Then, if `i` is not given it defaults to 0 for `k > 0` and `n - 1` for `k < 0`. If `j` is not given it defaults to `n` for `k > 0` and `-n-1` for `k < 0`. If `k` is not given it defaults to 1. Note that ``::`` is the same as ``:`` and means select all indices along this axis. From the above example:


.. code-block:: pycon

    >>> x[5:]
    array([5, 6, 7, 8, 9])

If the number of objects in the selection tuple is less than `N`, then ``:`` is assumed for any subsequent dimensions. For example:


.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm
    >>> x = np.array([[[1],[2],[3]], [[4],[5],[6]]])
    >>> x.shape
    (2, 3, 1)
    >>> x[1:2]
    array([[[4],
        [5],
        [6]]])

An integer, `i`, returns the same values as ``i:i+1`` except the dimensionality of the returned object is reduced by 1. In particular, a selection tuple with the p-th element an integer (and all other entries :) returns the corresponding sub-array with dimension `N - 1`. If `N = 1` then the returned object is an array scalar. 

If the selection tuple has all entries : except the p-th entry which is a slice object ``i:j:k``, then the returned array has dimension `N` formed by concatenating the sub-arrays returned by integer indexing of elements `i, i+k, …, i + (m - 1) k < j`.

Basic slicing with more than one non-``:`` entry in the slicing tuple, acts like repeated application of slicing using a single non-``:`` entry, where the non-``:`` entries are successively taken (with all other non-``:`` entries replaced by :). Thus, ``x[ind1, ..., ind2,:]`` acts like ``x[ind1][..., ind2, :]`` under basic slicing.

You may use slicing to set values in the array, but (unlike lists) you can never grow the array. The size of the value to be set in ``x[obj] = value`` must be (broadcastable to) the same shape as ``x[obj]``.

A slicing tuple can always be constructed as `obj` and used in the ``x[obj]`` notation. Slice objects can be used in the construction in place of the ``[start:stop:step]`` notation. For example, ``x[1:10:5, ::-1]`` can also be implemented as ``obj = (slice(1, 10, 5), slice(None, None, -1)); x[obj]``. This can be useful for constructing generic code that works on arrays of arbitrary dimensions.




Dimensional indexing tools
----------------------------------------------------------------


There are some tools to facilitate the easy matching of array shapes with expressions and in assignments.

Ellipsis expands to the number of ``:`` objects needed for the selection tuple to index all dimensions. In most cases, this means that the length of the expanded selection tuple is x.ndim. There may only be a single ellipsis present. From the above example:


.. code-block:: pycon

    >>> x[..., 0]
    array([[1, 2, 3],
      [4, 5, 6]])

This is equivalent to:

.. code-block:: pycon

    >>> x[:, :, 0]
    array([[1, 2, 3],
      [4, 5, 6]])

Each newaxis object in the selection tuple serves to expand the dimensions of the resulting selection by one unit-length dimension. The added dimension is the position of the newaxis object in the selection tuple. newaxis is an alias for None, and None can be used in place of this with the same result. From the above example:


.. code-block:: pycon

    >>> x[:, np.newaxis, :, :].shape
    (2, 1, 3, 1)
    >>> x[:, None, :, :].shape
    (2, 1, 3, 1)

This can be handy to combine two arrays in a way that otherwise would require explicit reshaping operations. For example:


.. code-block:: pycon

    >>> x = np.arange(5)
    >>> x[:, np.newaxis] + x[np.newaxis, :]
    array([[0, 1, 2, 3, 4],
      [1, 2, 3, 4, 5],
      [2, 3, 4, 5, 6],
      [3, 4, 5, 6, 7],
      [4, 5, 6, 7, 8]])









Broadcasting Overview
----------------------------------------------------------------


https://numpy.org/doc/stable/user/basics.broadcasting.html#basics-broadcasting

The term broadcasting describes how NumPy treats arrays with different shapes during arithmetic operations. Subject to certain constraints, the smaller array is “broadcast” across the larger array so that they have compatible shapes. Broadcasting provides a means of vectorizing array operations so that looping occurs in C instead of Python. It does this without making needless copies of data and usually leads to efficient algorithm implementations. There are, however, cases where broadcasting is a bad idea because it leads to inefficient use of memory that slows computation.

NumPy operations are usually done on pairs of arrays on an element-by-element basis. In the simplest case, the two arrays must have exactly the same shape, as in the following example:




.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm
    >>> a = np.array([1.0, 2.0, 3.0])
    >>> b = np.array([2.0, 2.0, 2.0])
    >>> a * b
    array([2.,  4.,  6.])

NumPy’s broadcasting rule relaxes this constraint when the arrays’ shapes meet certain constraints. The simplest broadcasting example occurs when an array and a scalar value are combined in an operation:


.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm
    >>> a = np.array([1.0, 2.0, 3.0])
    >>> b = 2.0
    >>> a * b
    array([2.,  4.,  6.])

The result is equivalent to the previous example where b was an array. We can think of the scalar b being stretched during the arithmetic operation into an array with the same shape as a. The new elements in b, as shown in Figure 1, are simply copies of the original scalar. The stretching analogy is only conceptual. NumPy is smart enough to use the original scalar value without actually making copies so that broadcasting operations are as memory and computationally efficient as possible.


FIGURE 1


The code in the second example is more efficient than that in the first because broadcasting moves less memory around during the multiplication (b is a scalar rather than an array).






General Broadcasting Rules
----------------------------------------------------------------


When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing (i.e. rightmost) dimension and works its way left. Two dimensions are compatible when they are equal, or one of them is 1.

If these conditions are not met, a ``ValueError: operands could not be broadcast together`` exception is thrown, indicating that the arrays have incompatible shapes.

Input arrays do not need to have the same *number* of dimensions. The resulting array will have the same number of dimensions as the input array with the greatest number of dimensions, where the size of each dimension is the largest size of the corresponding dimension among the input arrays. Note that missing dimensions are assumed to have size one.

For example, if you have a 256x256x3 array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:


.. code-block:: pycon

    Image  (3d array): 256 x 256 x 3
    Scale  (1d array):             3
    Result (3d array): 256 x 256 x 3

When either of the dimensions compared is one, the other is used. In other words, dimensions with size 1 are stretched or "copied" to match the other.

In the following example, both the A and B arrays have axes with length one that are expanded to a larger size during the broadcast operation:


.. code-block:: pycon

    A      (4d array):  8 x 1 x 6 x 1
    B      (3d array):      7 x 1 x 5
    Result (4d array):  8 x 7 x 6 x 5







Broadcastable arrays
----------------------------------------------------------------


A set of arrays is called "broadcastable" to the same shape if the above rules produce a valid result.

For example, if a.shape is (5,1), b.shape is (1,6), c.shape is (6,) and d.shape is () so that d is a scalar, then a, b, c, and d are all broadcastable to dimension (5,6); and

a acts like a (5,6) array where a[:,0] is broadcast to the other columns,

b acts like a (5,6) array where b[0,:] is broadcast to the other rows,

c acts like a (1,6) array and therefore like a (5,6) array where c[:] is broadcast to every row, and finally,

d acts like a (5,6) array where the single value is repeated.

Here are some more examples:


.. code-block:: pycon

    A      (2d array):  5 x 4
    B      (1d array):      1
    Result (2d array):  5 x 4

    A      (2d array):  5 x 4
    B      (1d array):      4
    Result (2d array):  5 x 4

    A      (3d array):  15 x 3 x 5
    B      (3d array):  15 x 1 x 5
    Result (3d array):  15 x 3 x 5

    A      (3d array):  15 x 3 x 5
    B      (2d array):       3 x 5
    Result (3d array):  15 x 3 x 5

    A      (3d array):  15 x 3 x 5
    B      (2d array):       3 x 1
    Result (3d array):  15 x 3 x 5


Here are examples of shapes that do not broadcast:

.. code-block:: python

    A      (1d array):  3
    B      (1d array):  4 # trailing dimensions do not match

    A      (2d array):      2 x 1
    B      (3d array):  8 x 4 x 3 # second from last dimensions mismatched


An example of broadcasting when a 1-d array is added to a 2-d array:


.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm

    >>> a = np.array([[ 0.0,  0.0,  0.0], [10.0, 10.0, 10.0], [20.0, 20.0, 20.0], [30.0, 30.0, 30.0]])
    >>> b = np.array([1.0, 2.0, 3.0])
    >>> a + b
    array([[  1.,   2.,   3.],
            [11.,  12.,  13.],
            [21.,  22.,  23.],
            [31.,  32.,  33.]])
    >>> b = np.array([1.0, 2.0, 3.0, 4.0])
    >>> a + b
    Traceback (most recent call last):
    ValueError: operands could not be broadcast together with shapes (4,3) (4,)



When the trailing dimensions of the arrays are unequal, broadcasting fails because it is impossible to align the values in the rows of the 1st array with the elements of the 2nd arrays for element-by-element addition.

Broadcasting provides a convenient way of taking the outer product (or any other outer operation) of two arrays. The following example shows an outer addition operation of two 1-d arrays:



.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm

    >>> a = np.array([0.0, 10.0, 20.0, 30.0])
    >>> b = np.array([1.0, 2.0, 3.0])
    >>> a[:, np.newaxis] + b
    array([[  1.,   2.,   3.],
            [11.,  12.,  13.],
            [21.,  22.,  23.],
            [31.,  32.,  33.]])



Here the newaxis index operator inserts a new axis into a, making it a two-dimensional 4x1 array. Combining the 4x1 array with b, which has shape (3,), yields a 4x3 array.









Assigning values to indexed arrays
----------------------------------------------------------------

Some text

https://numpy.org/doc/stable/user/basics.indexing.html#assigning-values-to-indexed-arrays

As mentioned, one can select a subset of an array to assign to using a single index, slices, and index and mask arrays. The value being assigned to the indexed array must be shape consistent (the same shape or broadcastable to the shape the index produces). For example, it is permitted to assign a constant to a slice:


.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm

    >>> x = np.arange(10)
    >>> x[2:7] = 1

or an array of the right size:


.. code-block:: pycon

    >>> x[2:7] = np.arange(5)

Note that assignments may result in changes if assigning higher types to lower types (like floats to ints) or even exceptions (assigning complex to floats or ints):



.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm

    >>> x[1] = 1.2
    >>> x[1]
    1
    >>> x[1] = 1.2j
    Traceback (most recent call last):
      ...
    TypeError: can't convert complex to int

Unlike some of the references (such as array and mask indices) assignments are always made to the original data in the array (indeed, nothing else would make sense!). Note though, that some actions may not work as one may naively expect. This particular example is often surprising to people:




.. code-block:: pycon

    >>> import numpy as np; from mpfunlab import mpm, dpm

    >>> x = np.arange(0, 50, 10)
    >>> x
    array([ 0, 10, 20, 30, 40])
    >>> x[np.array([1, 1, 3, 1])] += 1
    >>> x
    array([ 0, 11, 20, 31, 40])

Where people expect that the 1st location will be incremented by 3. In fact, it will only be incremented by 1. The reason is that a new array is extracted from the original (as a temporary) containing the values at 1, 1, 3, 1, then the value 1 is added to the temporary, and then the temporary is assigned back to the original array. Thus the value of the array at x[1] + 1 is assigned to x[1] three times, rather than being incremented 3 times.