<commit_before>.. currentmodule:: numpy

*************************
Numpy C Code Explanations
*************************

    Fanaticism consists of redoubling your efforts when you have forgotten
    your aim.
    --- *<NAME>*

    An authority is a person who can tell you more about something than
    you really care to know.
    --- *Unknown*

This Chapter attempts to explain the logic behind some of the new
pieces of code. The purpose behind these explanations is to enable
somebody to be able to understand the ideas behind the implementation
somewhat more easily than just staring at the code. Perhaps in this
way, the algorithms can be improved on, borrowed from, and/or
optimized.


Memory model
============

.. index::
   pair: ndarray; memory model

One fundamental aspect of the ndarray is that an array is seen as a
"chunk" of memory starting at some location. The interpretation of
this memory depends on the stride information. For each dimension in
an :math:`N` -dimensional array, an integer (stride) dictates how many
bytes must be skipped to get to the next element in that dimension.
Unless you have a single-segment array, this stride information must
be consulted when traversing through an array. It is not difficult to
write code that accepts strides, you just have to use (char \*)
pointers because strides are in units of bytes. Keep in mind also that
strides do not have to be unit-multiples of the element size. Also,
remember that if the number of dimensions of the array is 0 (sometimes
called a rank-0 array), then the strides and dimensions variables are
NULL.

Besides the structural information contained in the strides and
dimensions members of the :c:type:`PyArrayObject`, the flags contain
important information about how the data may be accessed. In particular,
the :c:data:`NPY_ARRAY_ALIGNED` flag is set when the memory is on a
suitable boundary according to the data-type array. Even if you have
a contiguous chunk of memory, you cannot just assume it is safe to
dereference a data- type-specific pointer to an element. Only if the
:c:data:`NPY_ARRAY_ALIGNED` flag is set is this a safe operation (on
some platforms it will work but on others, like Solaris, it will cause
a bus error). The :c:data:`NPY_ARRAY_WRITEABLE` should also be ensured
if you plan on writing to the memory area of the array. It is also
possible to obtain a pointer to an unwritable memory area. Sometimes,
writing to the memory area when the :c:data:`NPY_ARRAY_WRITEABLE` flag is not
set will just be rude. Other times it can cause program crashes ( *e.g.*
a data-area that is a read-only memory-mapped file).


Data-type encapsulation
=======================

.. index::
   single: dtype

The data-type is an important abstraction of the ndarray. Operations
will look to the data-type to provide the key functionality that is
needed to operate on the array. This functionality is provided in the
list of function pointers pointed to by the 'f' member of the
:c:type:`PyArray_Descr` structure. In this way, the number of data-types can be
extended simply by providing a :c:type:`PyArray_Descr` structure with suitable
function pointers in the 'f' member. For built-in types there are some
optimizations that by-pass this mechanism, but the point of the data-
type abstraction is to allow new data-types to be added.

One of the built-in data-types, the void data-type allows for
arbitrary structured types containing 1 or more fields as elements of the
array. A field is simply another data-type object along with an offset
into the current structured type. In order to support arbitrarily nested
fields, several recursive implementations of data-type access are
implemented for the void type. A common idiom is to cycle through the
elements of the dictionary and perform a specific operation based on
the data-type object stored at the given offset. These offsets can be
arbitrary numbers. Therefore, the possibility of encountering mis-
aligned data must be recognized and taken into account if necessary.


N-D Iterators
=============

.. index::
   single: array iterator

A very common operation in much of NumPy code is the need to iterate
over all the elements of a general, strided, N-dimensional array. This
operation of a general-purpose N-dimensional loop is abstracted in the
notion of an iterator object. To write an N-dimensional loop, you only
have to create an iterator object from an ndarray, work with the
dataptr member of the iterator object structure and call the macro
:c:func:`PyArray_ITER_NEXT` (it) on the iterator object to move to the next
element. The "next" element is always in C-contiguous order. The macro
works by first special casing the C-contiguous, 1-D, and 2-D cases
which work very simply.

For the general case, the iteration works by keeping track of a list
of coordinate counters in the iterator object. At each iteration, the
last coordinate counter is increased (starting from 0). If this
counter is smaller then one less than the size of the array in that
dimension (a pre-computed and stored value), then the counter is
increased and the dataptr member is increased by the strides in that
dimension and the macro ends. If the end of a dimension is reached,
the counter for the last dimension is reset to zero and the dataptr is
moved back to the beginning of that dimension by subtracting the
strides value times one less than the number of elements in that
dimension (this is also pre-computed and stored in the backstrides
member of the iterator object). In this case, the macro does not end,
but a local dimension counter is decremented so that the next-to-last
dimension replaces the role that the last dimension played and the
previously-described tests are executed again on the next-to-last
dimension. In this way, the dataptr is adjusted appropriately for
arbitrary striding.

The coordinates member of the :c:type:`PyArrayIterObject` structure maintains
the current N-d counter unless the underlying array is C-contiguous in
which case the coordinate counting is by-passed. The index member of
the :c:type:`PyArrayIterObject` keeps track of the current flat index of the
iterator. It is updated by the :c:func:`PyArray_ITER_NEXT` macro.


Broadcasting
============

.. index::
   single: broadcasting

In Numeric, broadcasting was implemented in several lines of code
buried deep in ufuncobject.c. In NumPy, the notion of broadcasting has
been abstracted so that it can be performed in multiple places.
Broadcasting is handled by the function :c:func:`PyArray_Broadcast`. This
function requires a :c:type:`PyArrayMultiIterObject` (or something that is a
binary equivalent) to be passed in. The :c:type:`PyArrayMultiIterObject` keeps
track of the broadcast number of dimensions and size in each
dimension along with the total size of the broadcast result. It also
keeps track of the number of arrays being broadcast and a pointer to
an iterator for each of the arrays being broadcast.

The :c:func:`PyArray_Broadcast` function takes the iterators that have already
been defined and uses them to determine the broadcast shape in each
dimension (to create the iterators at the same time that broadcasting
occurs then use the :c:func:`PyMultiIter_New` function). Then, the iterators are
adjusted so that each iterator thinks it is iterating over an array
with the broadcast size. This is done by adjusting the iterators
number of dimensions, and the shape in each dimension. This works
because the iterator strides are also adjusted. Broadcasting only
adjusts (or adds) length-1 dimensions. For these dimensions, the
strides variable is simply set to 0 so that the data-pointer for the
iterator over that array doesn't move as the broadcasting operation
operates over the extended dimension.

Broadcasting was always implemented in Numeric using 0-valued strides
for the extended dimensions. It is done in exactly the same way in
NumPy. The big difference is that now the array of strides is kept
track of in a :c:type:`PyArrayIterObject`, the iterators involved in a
broadcast result are kept track of in a :c:type:`PyArrayMultiIterObject`,
and the :c:func:`PyArray_BroadCast` call implements the broad-casting rules.


Array Scalars
=============

.. index::
   single: array scalars

The array scalars offer a hierarchy of Python types that allow a one-
to-one correspondence between the data-type stored in an array and the
Python-type that is returned when an element is extracted from the
array. An exception to this rule was made with object arrays. Object
arrays are heterogeneous collections of arbitrary Python objects. When
you select an item from an object array, you get back the original
Python object (and not an object array scalar which does exist but is
rarely used for practical purposes).

The array scalars also offer the same methods and attributes as arrays
with the intent that the same code can be used to support arbitrary
dimensions (including 0-dimensions). The array scalars are read-only
(immutable) with the exception of the void scalar which can also be
written to so that structured array field setting works more naturally
(a[0]['f1'] = ``value`` ).


Indexing
========

.. index::
   single: indexing

All python indexing operations ``arr[index]`` are organized by first preparing
the index and finding the index type. The supported index types are:
* integer
* newaxis
* slice
* ellipsis
* integer arrays/array-likes (fancy)
* boolean (single boolean array); if there is more than one boolean array as
  index or the shape does not match exactly, the boolean array will be
  converted to an integer array instead.
* 0-d boolean (and also integer); 0-d boolean arrays are a special
  case which has to be handled in the advanced indexing code. They signal
  that a 0-d boolean array had to be interpreted as an integer array.

As well as the scalar array special case signaling that an integer array
was interpreted as an integer index, which is important because an integer
array index forces a copy but is ignored if a scalar is returned (full integer
index). The prepared index is guaranteed to be valid with the exception of
out of bound values and broadcasting errors for advanced indexing. This
includes that an ellipsis is added for incomplete indices for example when
a two dimensional array is indexed with a single integer.

The next step depends on the type of index which was found. If all
dimensions are indexed with an integer a scalar is returned or set. A
single boolean indexing array will call specialized boolean functions.
Indices containing an ellipsis or slice but no advanced indexing will
always create a view into the old array by calculating the new strides and
memory offset.  This view can then either be returned or, for assignments,
filled using :c:func:`PyArray_CopyObject`. Note that `PyArray_CopyObject`
may also be called on temporary arrays in other branches to support
complicated assignments when the array is of object dtype.

Advanced indexing
-----------------

By far the most complex case is advanced indexing, which may or may not be
combined with typical view based indexing. Here integer indices are
interpreted as view based. Before trying to understand this, you may want
to make yourself familiar with its subtleties. The advanced indexing code
has three different branches and one special case:
* There is one indexing array and it, as well as the assignment array, can
  be iterated trivially. For example they may be contiguous. Also the
  indexing array must be of `intp` type and the value array in assignments
  should be of the correct type. This is purely a fast path.
* There are only integer array indices so that no subarray exists.
* View based and advanced indexing is mixed. In this case the view based
  indexing defines a collection of subarrays that are combined by the<commit_msg>DOC: Fix bad rst list formats in internals.code-explanations.rst.<commit_after>.. currentmodule:: numpy

*************************
Numpy C Code Explanations
*************************

    Fanaticism consists of redoubling your efforts when you have forgotten
    your aim.
    --- *<NAME>*

    An authority is a person who can tell you more about something than
    you really care to know.
    --- *Unknown*

This Chapter attempts to explain the logic behind some of the new
pieces of code. The purpose behind these explanations is to enable
somebody to be able to understand the ideas behind the implementation
somewhat more easily than just staring at the code. Perhaps in this
way, the algorithms can be improved on, borrowed from, and/or
optimized.


Memory model
============

.. index::
   pair: ndarray; memory model

One fundamental aspect of the ndarray is that an array is seen as a
"chunk" of memory starting at some location. The interpretation of
this memory depends on the stride information. For each dimension in
an :math:`N` -dimensional array, an integer (stride) dictates how many
bytes must be skipped to get to the next element in that dimension.
Unless you have a single-segment array, this stride information must
be consulted when traversing through an array. It is not difficult to
write code that accepts strides, you just have to use (char \*)
pointers because strides are in units of bytes. Keep in mind also that
strides do not have to be unit-multiples of the element size. Also,
remember that if the number of dimensions of the array is 0 (sometimes
called a rank-0 array), then the strides and dimensions variables are
NULL.

Besides the structural information contained in the strides and
dimensions members of the :c:type:`PyArrayObject`, the flags contain
important information about how the data may be accessed. In particular,
the :c:data:`NPY_ARRAY_ALIGNED` flag is set when the memory is on a
suitable boundary according to the data-type array. Even if you have
a contiguous chunk of memory, you cannot just assume it is safe to
dereference a data- type-specific pointer to an element. Only if the
:c:data:`NPY_ARRAY_ALIGNED` flag is set is this a safe operation (on
some platforms it will work but on others, like Solaris, it will cause
a bus error). The :c:data:`NPY_ARRAY_WRITEABLE` should also be ensured
if you plan on writing to the memory area of the array. It is also
possible to obtain a pointer to an unwritable memory area. Sometimes,
writing to the memory area when the :c:data:`NPY_ARRAY_WRITEABLE` flag is not
set will just be rude. Other times it can cause program crashes ( *e.g.*
a data-area that is a read-only memory-mapped file).


Data-type encapsulation
=======================

.. index::
   single: dtype

The data-type is an important abstraction of the ndarray. Operations
will look to the data-type to provide the key functionality that is
needed to operate on the array. This functionality is provided in the
list of function pointers pointed to by the 'f' member of the
:c:type:`PyArray_Descr` structure. In this way, the number of data-types can be
extended simply by providing a :c:type:`PyArray_Descr` structure with suitable
function pointers in the 'f' member. For built-in types there are some
optimizations that by-pass this mechanism, but the point of the data-
type abstraction is to allow new data-types to be added.

One of the built-in data-types, the void data-type allows for
arbitrary structured types containing 1 or more fields as elements of the
array. A field is simply another data-type object along with an offset
into the current structured type. In order to support arbitrarily nested
fields, several recursive implementations of data-type access are
implemented for the void type. A common idiom is to cycle through the
elements of the dictionary and perform a specific operation based on
the data-type object stored at the given offset. These offsets can be
arbitrary numbers. Therefore, the possibility of encountering mis-
aligned data must be recognized and taken into account if necessary.


N-D Iterators
=============

.. index::
   single: array iterator

A very common operation in much of NumPy code is the need to iterate
over all the elements of a general, strided, N-dimensional array. This
operation of a general-purpose N-dimensional loop is abstracted in the
notion of an iterator object. To write an N-dimensional loop, you only
have to create an iterator object from an ndarray, work with the
dataptr member of the iterator object structure and call the macro
:c:func:`PyArray_ITER_NEXT` (it) on the iterator object to move to the next
element. The "next" element is always in C-contiguous order. The macro
works by first special casing the C-contiguous, 1-D, and 2-D cases
which work very simply.

For the general case, the iteration works by keeping track of a list
of coordinate counters in the iterator object. At each iteration, the
last coordinate counter is increased (starting from 0). If this
counter is smaller then one less than the size of the array in that
dimension (a pre-computed and stored value), then the counter is
increased and the dataptr member is increased by the strides in that
dimension and the macro ends. If the end of a dimension is reached,
the counter for the last dimension is reset to zero and the dataptr is
moved back to the beginning of that dimension by subtracting the
strides value times one less than the number of elements in that
dimension (this is also pre-computed and stored in the backstrides
member of the iterator object). In this case, the macro does not end,
but a local dimension counter is decremented so that the next-to-last
dimension replaces the role that the last dimension played and the
previously-described tests are executed again on the next-to-last
dimension. In this way, the dataptr is adjusted appropriately for
arbitrary striding.

The coordinates member of the :c:type:`PyArrayIterObject` structure maintains
the current N-d counter unless the underlying array is C-contiguous in
which case the coordinate counting is by-passed. The index member of
the :c:type:`PyArrayIterObject` keeps track of the current flat index of the
iterator. It is updated by the :c:func:`PyArray_ITER_NEXT` macro.


Broadcasting
============

.. index::
   single: broadcasting

In Numeric, broadcasting was implemented in several lines of code
buried deep in ufuncobject.c. In NumPy, the notion of broadcasting has
been abstracted so that it can be performed in multiple places.
Broadcasting is handled by the function :c:func:`PyArray_Broadcast`. This
function requires a :c:type:`PyArrayMultiIterObject` (or something that is a
binary equivalent) to be passed in. The :c:type:`PyArrayMultiIterObject` keeps
track of the broadcast number of dimensions and size in each
dimension along with the total size of the broadcast result. It also
keeps track of the number of arrays being broadcast and a pointer to
an iterator for each of the arrays being broadcast.

The :c:func:`PyArray_Broadcast` function takes the iterators that have already
been defined and uses them to determine the broadcast shape in each
dimension (to create the iterators at the same time that broadcasting
occurs then use the :c:func:`PyMultiIter_New` function). Then, the iterators are
adjusted so that each iterator thinks it is iterating over an array
with the broadcast size. This is done by adjusting the iterators
number of dimensions, and the shape in each dimension. This works
because the iterator strides are also adjusted. Broadcasting only
adjusts (or adds) length-1 dimensions. For these dimensions, the
strides variable is simply set to 0 so that the data-pointer for the
iterator over that array doesn't move as the broadcasting operation
operates over the extended dimension.

Broadcasting was always implemented in Numeric using 0-valued strides
for the extended dimensions. It is done in exactly the same way in
NumPy. The big difference is that now the array of strides is kept
track of in a :c:type:`PyArrayIterObject`, the iterators involved in a
broadcast result are kept track of in a :c:type:`PyArrayMultiIterObject`,
and the :c:func:`PyArray_BroadCast` call implements the broad-casting rules.


Array Scalars
=============

.. index::
   single: array scalars

The array scalars offer a hierarchy of Python types that allow a one-
to-one correspondence between the data-type stored in an array and the
Python-type that is returned when an element is extracted from the
array. An exception to this rule was made with object arrays. Object
arrays are heterogeneous collections of arbitrary Python objects. When
you select an item from an object array, you get back the original
Python object (and not an object array scalar which does exist but is
rarely used for practical purposes).

The array scalars also offer the same methods and attributes as arrays
with the intent that the same code can be used to support arbitrary
dimensions (including 0-dimensions). The array scalars are read-only
(immutable) with the exception of the void scalar which can also be
written to so that structured array field setting works more naturally
(a[0]['f1'] = ``value`` ).


Indexing
========

.. index::
   single: indexing

All python indexing operations ``arr[index]`` are organized by first preparing
the index and finding the index type. The supported index types are:

* integer
* newaxis
* slice
* ellipsis
* integer arrays/array-likes (fancy)
* boolean (single boolean array); if there is more than one boolean array as
  index or the shape does not match exactly, the boolean array will be
  converted to an integer array instead.
* 0-d boolean (and also integer); 0-d boolean arrays are a special
  case which has to be handled in the advanced indexing code. They signal
  that a 0-d boolean array had to be interpreted as an integer array.

As well as the scalar array special case signaling that an integer array
was interpreted as an integer index, which is important because an integer
array index forces a copy but is ignored if a scalar is returned (full integer
index). The prepared index is guaranteed to be valid with the exception of
out of bound values and broadcasting errors for advanced indexing. This
includes that an ellipsis is added for incomplete indices for example when
a two dimensional array is indexed with a single integer.

The next step depends on the type of index which was found. If all
dimensions are indexed with an integer a scalar is returned or set. A
single boolean indexing array will call specialized boolean functions.
Indices containing an ellipsis or slice but no advanced indexing will
always create a view into the old array by calculating the new strides and
memory offset.  This view can then either be returned or, for assignments,
filled using :c:func:`PyArray_CopyObject`. Note that `PyArray_CopyObject`
may also be called on temporary arrays in other branches to support
complicated assignments when the array is of object dtype.

Advanced indexing
-----------------

By far the most complex case is advanced indexing, which may or may not be
combined with typical view based indexing. Here integer indices are
interpreted as view based. Before trying to understand this, you may want
to make yourself familiar with its subtleties. The advanced indexing code
has three different branches and one special case:

* There is one indexing array and it, as well as the assignment array, can
  be iterated trivially. For example they may be contiguous. Also the
  indexing array must be of `intp` type and the value array in assignments
  should be of the correct type. This is purely a fast path.
* There are only integer array indices so that no subarray exists.
* View based and advanced indexing is mixed. In this case the view based
  indexing defines a collection of subarrays that are combined by the
{'fraction_non_alphanumeric': 0.049806332098349615, 'fraction_numerical': 0.001263051532502526, 'mean_word_length': 4.488216266173752, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1687059', 'n_tokens_mistral': 5877, 'n_tokens_neox': 5606, 'n_words': 3688}