# 06 - Matrix, vector, cube and field classes

Adapted from the official Armadillo documentation.

# 1 Dense matrix class

`Mat<type>`, `mat` and `cx_mat` are classes for dense matrices, with elements stored in column-major ordering (e.g., column by column).

The root matrix class is `Mat<type>`, where `type` is one of:

- `float`
- `double`
- `std::complex<float>`
- `std::complex<double>`
- `short`
- `int`
- `long`
- `unsigned short`
- `unsigned int`
- `unsigned long`

For convenience the following typedefs have been defined:

- `mat = Mat<double>`
- `dmat = Mat<double>`
- `fmat = Mat<float>`
- `cx_mat = Mat<cx_double>` (`cx_double` is a shortcut for `std::complex<double>`)
- `cx_dmat = Mat<cx_double>`
- `cx_fmat = Mat<cx_float>` (`cx_float` is a shortcut for `std::complex<float>`)
- `umat = Mat<uword>` (`uword` is a shortcut for `unsigned int`)
- `imat = Mat<sword>` (`sword` is a shortcut for `signed int`)

The `mat` type is used for convenience, and it is possible to use other matrix types (e.g, `fmat`, `cx_mat`) instead.

Matrix types with integer elements (such as `umat` and `imat`) cannot hold special values such as NaN and Inf.

Functions which use LAPACK (generally matrix decompositions) are only valid for the following matrix types: `mat`, `dmat`, `fmat`, `cx_mat`, `cx_dmat`, `cx_fmat`.

## 1.1 Constructors

- `mat()`
- `mat(n_rows, n_cols)`
- `mat(n_rows, n_cols, fill_form)` (elements are initialised according to `fill_form`)
- `mat(size(X))`
- `mat(size(X), fill_form)` (elements are initialised according to `fill_form`)
- `mat(mat)`
- `mat(vec)`
- `mat(rowvec)`
- `mat(initializer_list)`
- `mat(string)`
- `mat(std::vector)` (treated as a column vector)
- `mat(sp_mat)` (for converting a sparse matrix to a dense matrix)
- `cx_mat(mat,mat)` (for constructing a complex matrix out of two real matrices)

The elements can be explicitly initialised during construction by specifying `fill_form`, which is one of:

- `fill::zeros` set all elements to 0 (default in armadillo4r)
- `fill::ones` set all elements to 1
- `fill::eye` set the elements on the main diagonal to 1 and off-diagonal elements to 0
- `fill::randu` set all elements to random values from a uniform distribution in the [0,1] interval
- `fill::randn` set all elements to random values from a normal distribution with zero mean and unit variance
- `fill::value(scalar)` set all elements to specified scalar
- `fill::none` do not initialise the elements (matrix may have garbage values)

For the `mat(string)` constructor, the format is elements separated by spaces, and rows denoted by semicolons. For example, the 2x2 identity matrix can be created using `"1 0; 0 1"`. Note that string based initialisation is slower than directly setting the elements or using element initialisation.

Each instance of `mat` automatically allocates and releases internal memory. All internally allocated memory used by an instance of mat is automatically released as soon as the instance goes out of scope. For example, if an instance of mat is declared inside a function, it will be automatically destroyed at the end of the function. To forcefully release memory at any point, use `.reset()`. Note that in normal use this is not required.

## 1.2 Advanced constructors

- `mat(ptr_aux_mem, n_rows, n_cols, copy_aux_mem = true, strict = false)`
  - Create a matrix using data from writable auxiliary (external) memory, where `ptr_aux_mem` is a pointer to the memory. By default the matrix allocates its own memory and copies data from the auxiliary memory (for safety). However, if `copy_aux_mem` is set to `false`, the matrix will instead directly use the auxiliary memory (e.g., no copying). This is faster, but can be dangerous unless you know what you are doing.
  - The `strict` parameter comes into effect only when copy_aux_mem is set to `false` (e.g., the matrix is directly using auxiliary memory).
    - When `strict` is set to `false`, the matrix will use the auxiliary memory until a size change or an aliasing event.
    - When `strict` is set to `true`, the matrix will be bound to the auxiliary memory for its lifetime. The number of elements in the matrix cannot be changed.
- `mat(const ptr_aux_mem, n_rows, n_cols)`
  - Create a matrix by copying data from read-only auxiliary memory, where `ptr_aux_mem` is a pointer to the memory
- `mat::fixed<n_rows, n_cols>`
  - Create a fixed size matrix, with the size specified via template arguments. Memory for the matrix is reserved at compile time. This is generally faster than dynamic memory allocation, but the size of the matrix cannot be changed afterwards (directly or indirectly).
  - For convenience, there are several pre-defined typedefs for each matrix type (where the types are: `umat`, `imat`, `fmat`, `mat`, `cx_fmat`, `cx_mat`). The typedefs specify a square matrix size, ranging from 2x2 to 9x9. The typedefs were defined by appending a two digit form of the size to the matrix type. Examples: `mat33` is equivalent to `mat::fixed<3,3>`, and `cx_mat44` is equivalent to `cx_mat::fixed<4,4>`.
- `mat::fixed<n_rows, n_cols>(fill_form)`
  - Create a fixed size matrix, with the elements explicitly initialised according to `fill_form`.
- `mat::fixed<n_rows, n_cols>(const ptr_aux_mem)`
  - Create a fixed size matrix, with the size specified via template arguments; data is copied from auxiliary memory, where ptr_aux_mem is a pointer to the memory.

## 1.3 Examples

    set.seed(123)
    a <- matrix(runif(25), nrow = 5, ncol = 5)

    [[cpp4r::register]] doubles_matrix<> matrix1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a); // convert from R to C++

      double x = A(0, 0);  // access an element on row 1, column 1
      A = A + x; // scalar addition

      mat B = A + A; // matrix addition
      mat C = A * B; // matrix multiplication
      mat D = A % B; // element-wise matrix multiplication

      mat res = B + C + D;

      return as_doubles_matrix(res);  // convert from C++ to R
    }

    [[cpp4r::register]] list matrix2_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);
      mat B = A + A;

      cx_mat X(A,B); // construct a complex matrix out of two real matrices

      B.zeros();                       // set all elements to zero
      B.set_size(A.n_rows, A.n_cols);  // resize the matrix
      B.ones(5, 6);                    // same as mat B(5, 6, fill::ones)

      mat::fixed<5,6> F; // fixed size matrix

      double aux_mem[24]; // auxiliary memory
      mat H(&aux_mem[0], 4, 6, false);  // use auxiliary memory

      X = X + F.submat(0, 0, 4, 4) + H(1, 2)

      Mat<double> res_real = real(X);
      Mat<double> res_imag = imag(X);

      writable::list res;
      res.push_back({"real"_nm = as_doubles_matrix(res_real)});
      res.push_back({"imag"_nm = as_doubles_matrix(res_imag)});

      return res;
    }

# 2 Column vector class

`Col<type>`, `vec` and `cx_vec` are classes for column vectors (dense matrices with one column).

The `Col<type>` class is derived from the `Mat<type>` class and inherits most of the member functions.

For convenience the following typedefs have been defined:

- `vec = colvec = Col<double>`
- `dvec = dcolvec = Col<double>`
- `fvec = fcolvec = Col<float>`
- `cx_vec = cx_colvec = Col<[cx_double](#cx_double)>`
- `cx_dvec = cx_dcolvec = Col<[cx_double](#cx_double)>`
- `cx_fvec = cx_fcolvec = Col<[cx_float](#cx_double)>`
- `uvec = ucolvec = Col<[uword](#uword)>`
- `ivec = icolvec = Col<[sword](#uword)>`

The `vec` and `colvec` types have the same meaning and are used interchangeably.

The types `vec` or `colvec` are used for convenience. It is possible to use other column vector types instead (e.g., `fvec` or `fcolvec`).

Functions which take `mat` as input can generally also take `Col` as input. Main exceptions are functions that require square matrices.

## 2.1 Constructors

- `vec()`
- `vec(_n_elem_)`
- `vec(_n_elem, fill_form)` (elements are initialised according to `fill_form`)
- `vec(size(X))`
- `vec(size(X), fill_form)` (elements are initialised according to `fill_form`)
- `vec(vec)`
- `vec(mat)` (`std::logic_error` exception is thrown if the given matrix has more than one column)
- `vec(initializer_list)`
- `vec(string)` (elements separated by spaces)
- `vec(std::vector)`
- `cx_vec(vec,vec)` (for constructing a complex vector out of two real vectors)

## 2.2 Advanced constructors

- `vec(ptr_aux_mem, number_of_elements, copy_aux_mem = true, strict = false)`:
  - Create a column vector using data from writable auxiliary (external) memory, where ptr_aux_mem is a pointer to the memory. By default the vector allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to `false`, the vector will instead directly use the auxiliary memory (e.g., no copying). This is faster, but can be dangerous unless you know what you are doing.
  - The `strict` parameter comes into effect only when `copy_aux_mem` is set to `false` (e.g., the vector is directly using auxiliary memory).
  - When `strict` is set to `false`, the vector will use the auxiliary memory until a size change or an aliasing event.
  - When `strict` is set to `true`, the vector will be bound to the auxiliary memory for its lifetime. The number of elements in the vector cannot be changed.
- `vec(const ptr_aux_mem, number_of_elements)`: Create a column vector by copying data from read-only auxiliary memory, where ptr_aux_mem is a pointer to the memory.
- `vec::fixed<number_of_elements>`:
  - Create a fixed size column vector, with the size specified via the template argument. Memory for the vector is reserved at compile time. This is generally faster than dynamic memory allocation, but the size of the vector cannot be changed afterwards (directly or indirectly).
  - For convenience, there are several pre-defined typedefs for each vector type (where the types are: `uvec`, `ivec`, `fvec`, `vec`, `cx_fvec`, `cx_vec` as well as the corresponding `colvec` versions). The pre-defined typedefs specify vector sizes ranging from 2 to 9. The typedefs were defined by appending a single digit form of the size to the vector type. Examples: `vec3` is equivalent to `vec::fixed<3>`, and `cx_vec4` is equivalent to `cx_vec::fixed<4>`.
- `vec::fixed<number_of_elements>(fill_form)`: Create a fixed size column vector, with the elements explicitly initialised according to `fill_form`.
- `vec::fixed<number_of_elements>(const ptr_aux_mem)`: Create a fixed size column vector, with the size specified via the template argument. The data is copied from auxiliary memory, where `ptr_aux_mem` is a pointer to the memory.

## 2.3 Examples

    set.seed(123)
    x <- runif(10)
    y <- rep(1, 10)

    [[cpp4r::register]] doubles column1_(const doubles& x, const doubles& y) {
      vec X = as_Col(x); // convert from R to C++
      vec Y = as_Col(y);

      mat A(10, 10, fill::randu);
      vec Z = A.col(5); // extract a column vector

      Z = Z + Y + X;

      return as_doubles(Z); // convert from C++ to R
    }

# 3 Row vector class

`Row<type>`, `rowvec` and `cx_rowvec` are classes for row vectors (dense matrices with one row).

The template `Row<type>` class is derived from the `Mat<type>` class and inherits most of the member functions.

For convenience the following typedefs have been defined:

- `rowvec = Row<double>`
- `drowvec = Row<double>`
- `frowvec = Row<float>`
- `cx_rowvec = Row<cx_double>`
- `cx_drowvec = Row<cx_double>`
- `cx_frowvec = Row<cx_float>`
- `urowvec = Row<uword>`
- `irowvec = Row<sword>`

The `rowvec` type is used for convenience. It is possible to use other row vector types instead (e.g., `frowvec`).

Functions which take `mat` as input can generally also take `Row` as input. Main exceptions are functions which require square matrices.

## 3.1 Constructors

- `rowvec()`
- `rowvec(n_elem)`
- `rowvec(n_elem, fill_form)` (elements are initialised according to `fill_form`)
- `rowvec(size(X))`
- `rowvec(size(X), fill_form)` (elements are initialised according to `fill_form`)
- `rowvec(rowvec)`
- `rowvec(mat)` (`std::logic_error` exception is thrown if the given matrix has more than one row)
- `rowvec(initializer_list)`
- `rowvec(string)` (elements separated by spaces)
- `rowvec(std::vector)`
- `cx_rowvec(rowvec,rowvec)` (for constructing a complex row vector out of two real row vectors)

## 3.2 Advanced constructors

- `rowvec(ptr_aux_mem, number_of_elements, copy_aux_mem = true, strict = false)`
  - Create a row vector using data from writable auxiliary (external) memory, where ptr_aux_mem is a pointer to the memory. By default the vector allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to `false`, the vector will instead directly use the auxiliary memory (e.g., no copying); this is faster, but can be dangerous unless you know what you are doing.
  - The `strict` parameter comes into effect only when copy_aux_mem is set to `false` (e.g., the vector is directly using auxiliary memory).
    - When `strict` is set to `false`, the vector will use the auxiliary memory until a size change or an aliasing event.
    - When `strict` is set to `true`, the vector will be bound to the auxiliary memory for its lifetime. The number of elements in the vector cannot be changed.
- `rowvec(const ptr_aux_mem, number_of_elements)`
  - Create a row vector by copying data from read-only auxiliary memory, where `ptr_aux_mem` is a pointer to the memory.
- `rowvec::fixed<number_of_elements>`
  - Create a fixed size row vector, with the size specified via the template argument. Memory for the vector is reserved at compile time. This is generally faster than dynamic memory allocation, but the size of the vector cannot be changed afterwards (directly or indirectly).
  - For convenience, there are several pre-defined typedefs for each vector type (where the types are: `urowvec`, `irowvec`, `frowvec`, `rowvec`, `cx_frowvec`, `cx_rowvec`). The pre-defined typedefs specify vector sizes ranging from 2 to 9. The typedefs were defined by appending a single digit form of the size to the vector type. Examples: `rowvec3` is equivalent to `rowvec::fixed<3>`, and `cx_rowvec4` is equivalent to `cx_rowvec::fixed<4>`.
- `rowvec::fixed<number_of_elements>(fill_form)`
  - Create a fixed size row vector, with the elements explicitly initialised according to `fill_form`.
- `rowvec::fixed<number_of_elements>(const ptr_aux_mem)`
  - Create a fixed size row vector, with the size specified via the template argument. The data is copied from auxiliary memory, where `ptr_aux_mem` is a pointer to the memory.

## 3.3 Examples

⚠️Important⚠️: ‘armadillo4r’ is an opinionated package and it follows the notation from Econometrics by Bruce E. Hansen. It intentionally exports/imports matrices and column vectors. You can use row vectors in the functions, but the communication between R and C++ does not accept row vectors unless you transpose or convert those to matrices.

    set.seed(123)
    x <- runif(10)
    y <- rep(1, 10)

    [[cpp4r::register]] doubles row1_(const doubles& x, const doubles& y) {
      vec X = as_Col(x);  // convert from R to C++
      vec Y = as_Col(y);

      mat A(10, 10, fill::randu);

      rowvec Z = A.row(5);  // extract a row vector
      Z = Z + Y.t() + X.t(); // transpose Y and X to be able to sum

      vec res = Z.t();

      return as_doubles(res);  // convert from C++ to R
    }

# 4 Cube class

`Cube<type>`, `cube` and `cx_cube` are classes for cubes, also known as quasi 3rd order tensors or “3D matrices”.

The data is stored as a set of slices (matrices) stored contiguously within memory. Within each slice, elements are stored with column-major ordering (e.g., column by column)

The root cube class is `Cube<type>`, where `type` is one of: `float`, `double`, `std::complex<float>`, `std::complex<double>`, `short`, `int`, `long` and unsigned versions of `short`, `int`, `long`.

For convenience the following typedefs have been defined:

- `cube = Cube<double>`
- `dcube = Cube<double>`
- `fcube = Cube<float>`
- `cx_cube = Cube<cx_double>`
- `cx_dcube = Cube<cx_double>`
- `cx_fcube = Cube<cx_float>`
- `ucube = Cube<uword>`
- `icube = Cube<sword>`

The `cube` type is used for convenience. It is possible to use other types instead (e.g., `fcube`).

Each cube slice can be interpreted as a matrix, hence functions which take `Mat` as input can generally also take cube slices as input.

## 4.1 Constructors

- `cube()`
- `cube(n_rows, n_cols, n_slices_)`
- `cube(n_rows, n_cols, n_slices, fill_form)` (elements are initialised according to `fill_form`)
- `cube(size(X))`
- `cube(size(X), fill_form)` (elements are initialised according to `fill_form`)
- `cube(cube)`
- `cx_cube(cube, cube)` (for constructing a complex cube out of two real cubes)

The elements can be explicitly initialised during construction by specifying `fill_form`, which is one of:

- `fill::zeros`: set all elements to 0 (default in armadillo4r)
- `fill::ones`: set all elements to 1
- `fill::randu`: set all elements to random values from a uniform distribution in the [0,1] interval
- `fill::randn`: set all elements to random values from a normal distribution with zero mean and unit variance
- `fill::value(scalar)`: set all elements to specified scalar
- `fill::none`: do not initialise the elements (cube may have garbage values)

Each instance of `cube` automatically allocates and releases internal memory. All internally allocated memory used by an instance of `cube` is automatically released as soon as the instance goes out of scope. For example, if an instance of `cube` is declared inside a function, it will be automatically destroyed at the end of the function. To forcefully release memory at any point, use `.reset()` note that in normal use this is not required.

## 4.2 Advanced constructors

- `cube::fixed<n_rows, n_cols, n_slices>`: Create a fixed size cube, with the size specified via template arguments. Memory for the cube is reserved at compile time. This is generally faster than dynamic memory allocation, but the size of the cube cannot be changed afterwards (directly or indirectly).
- `cube(ptr_aux_mem, n_rows, n_cols, n_slices, copy_aux_mem = true, strict = false)`:
  - Create a cube using data from writable auxiliary (external) memory, where `ptr_aux_mem` is a pointer to the memory. By default the cube allocates its own memory and copies data from the auxiliary memory (for safety). However, if `copy_aux_mem` is set to `false`, the cube will instead directly use the auxiliary memory (e.g., no copying). This is faster, but can be dangerous unless you know what you are doing.
  - The `strict` parameter comes into effect only when `copy_aux_mem` is set to `false` (e.g., the cube is directly using auxiliary memory).
  - When `strict` is set to `false`, the cube will use the auxiliary memory until a size change or an aliasing event.
  - When `strict` is set to `true`, the cube will be bound to the auxiliary memory for its lifetime. The number of elements in the cube cannot be changed.
- `cube(const ptr_aux_mem, n_rows, n_cols, n_slices)`: Create a cube by copying data from read-only auxiliary memory, where `ptr_aux_mem` is a pointer to the memory.

## 4.3 Examples

    set.seed(123)
    a <- matrix(runif(4), nrow = 2, ncol = 2)
    b <- matrix(rnorm(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> cube1_(const doubles_matrix<>& a,
                                                    const doubles_matrix<>& b) {
      mat A = as_Mat(a);  // convert from R to C++
      mat B = as_Mat(b);

      cube X(A.n_rows, A.n_cols, 2);  // create a cube with 2 slices
      X.slice(0) = A;                 // copy A into first slice
      X.slice(1) = B;                 // copy B into second slice

      cube Y = X + X;  // cube addition
      cube Z = X % X;  // element-wise cube multiplication

      mat res = Y.slice(0) + Z.slice(1);

      return as_doubles_matrix(res);  // convert from C++ to R
    }

The size of individual slices cannot be changed. The following will not work:

    cube c(5,6,7);
    c.slice(0) = randu<mat>(10,20); // wrong size

# 5 Field class

`field<object_type>` is a class for storing arbitrary objects in matrix-like or cube-like layouts.

It is similar to a matrix or cube, but instead of each element being a scalar, each element can be a vector, or matrix, or cube. This is similar to a list in R.

Each element can have an arbitrary size (e.g., in a field of matrices, each matrix can have a unique size).

## 5.1 Constructors

`object_type` is another class (e.g., `vec`, `mat`, `std::string`, etc)

- `field<object_type>()`
- `field<object_type>(n_elem)`
- `field<object_type>(n_rows, n_cols)`
- `field<object_type>(n_rows, n_cols, n_slices)`
- `field<object_type>(size(X))`
- `field<object_type>(field<object_type>)`

## 5.2 Examples

    set.seed(123)
    a <- matrix(runif(4), nrow = 2, ncol = 2)
    b <- matrix(rnorm(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> field1_(const doubles_matrix<>& a,
                                                     const doubles_matrix<>& b) {
      mat A = as_Mat(a);  // convert from R to C++
      mat B = as_Mat(b);

      field<mat> F(A.n_rows, A.n_cols, 3);  // create a field with 2 matrices
      F(0) = A;                             // copy A into first location
      F(1) = B;                             // copy B into second location
      F(2) = F(0) + F(1);                   // matrix addition

      mat res = F(0) + F(1) + F(2).t();

      return as_doubles_matrix(res);  // convert from C++ to R
    }

## 5.3 Caveat

To store a set of matrices of the same size, the `Cube` class is more efficient.

# 6 Member Functions and Variables

Function/Variable | Description
--- | ---
`.n_rows` | number of rows
`.n_cols` | number of columns
`.n_elem` | number of elements
`.n_slices` | number of slices
`()` | element access
`[]` | element access
`.at()` | element access
`.zeros` | set all elements to zero
`.ones` | set all elements to one
`.eye` | set elements along main diagonal to one and off-diagonal elements to zero
`.randu` | set all elements to random values from a uniform distribution
`.randn` | set all elements to random values from a normal distribution
`.fill` | set all elements to specified value
`.imbue` | imbue (fill) with values provided by functor or lambda function
`.clean` | replace elements below a threshold with zeros
`.replace` | replace specific elements with a new value
`.clamp` | clamp values to lower and upper limits
`.transform` | transform each element via functor or lambda function
`.for_each` | apply a functor or lambda function to each element
`.set_size` | change size without keeping elements (fast)
`.reshape` | change size while keeping elements
`.resize` | change size while keeping elements and preserving layout
`.copy_size` | change size to be same as given object
`.reset` | change size to empty
`.diag` | read/write access to matrix diagonals
`.each_col` | vector operations applied to each column of matrix (aka “broadcasting”)
`.each_row` | vector operations applied to each row of matrix (aka “broadcasting”)
`.each_slice` | matrix operations applied to each slice of cube (aka “broadcasting”)
`.set_imag` | set imaginary part
`.set_real` | set real part
`.insert_rows` | insert vector/matrix/cube at specified row
`.insert_cols` | insert vector/matrix/cube at specified column
`.insert_slices` | insert vector/matrix/cube at specified slice
`.shed_rows` | remove specified rows
`.shed_cols` | remove specified columns
`.shed_slices` | remove specified slices
`.swap_rows` | swap specified rows
`.swap_cols` | swap specified columns
`.swap` | swap contents with given object
`.memptr` | raw pointer to memory
`.colptr` | raw pointer to memory used by specified column
`.as_col` | return flattened matrix column as column vector
`.as_row` | return flattened matrix row as row vector
`.col_as_mat` | return matrix representation of cube column
`.row_as_mat` | return matrix representation of cube row
`.as_dense` | return dense vector/matrix representation of sparse matrix expression
`.t` | return matrix transpose
`.st` | return matrix conjugate transpose
`.i` | return inverse of square matrix
`.min` | return minimum value
`.max` | return maximum value
`.index_min` | return index of minimum value
`.index_max` | return index of maximum value
`.eval` | force evaluation of delayed expression
`.in_range` | check whether given location or span is valid
`.is_empty` | check whether object is empty
`.is_vec` | check whether matrix is a vector
`.is_sorted` | check whether vector or matrix is sorted
`.is_trimatu` | check whether matrix is upper triangular
`.is_trimatl` | check whether matrix is lower triangular
`.is_diagmat` | check whether matrix is diagonal
`.is_square` | check whether matrix is square sized
`.is_symmetric` | check whether matrix is symmetric
`.is_hermitian` | check whether matrix is hermitian
`.is_sympd` | check whether matrix is symmetric/hermitian positive definite
`.is_zero` | check whether all elements are zero
`.is_finite` | check whether all elements are finite
`.has_inf` | check whether any element is +/- infinity
`.has_nan` | check whether any element is NaN

# 7 Attributes

`n_*` provides information for different objects:

- `.n_rows` number of rows for `Mat`, `Col`, `Row`, `Cube`, `field`, and `SpMat`.
- `.n_cols` number of columns for `Mat`, `Col`, `Row`, `Cube`, `field`, and `SpMat`.
- `.n_elem` total number of elements for `Mat`, `Col`, `Row`, `Cube`, `field`, and `SpMat`.
- `.n_slices` number of slices for `Cube` and `field`.
- `.n_nonzero` number of non-zero elements for `SpMat`.

For the `Col` and `Row` classes, `n_elem` also indicates vector length.

The variables are read-only and of type `uword`. To change the size, use `set_size`, `copy_size`, `zeros_member`, `ones_member`, or `reset`.

To avoid compiler warnings about implicit conversion when operating `uword` with `integers`/`doubles` to pass data to R, converte `uword` to `int` with `static_cast<int>` or declare these as `int`.

## 7.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] integers attr1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++

      // uword or int can be used
      int n_rows = A.n_rows;  // number of rows
      int n_cols = A.n_cols;  // number of columns
      int n_elem = A.n_elem;  // number of elements

      writable::integers res({n_rows, n_cols, n_elem});
      res.attr("names") = strings({"n_rows", "n_cols", "n_elem"});

      return res;
    }

# 8 Element/object access

Provide access to individual elements or objects stored in a container object (e.g., `Mat`, `Col`, `Row`, `Cube`, `field`).

- `(i)` For `vec` and `rowvec`, access the element stored at index `i`. For `Mat`, `Cube` and `field`, access the element/object stored at index `i` under the assumption of a flat layout, with column-major ordering of data (e.g., column by column). An exception is thrown if the requested element is out of bounds.
- `.at(i)` or `[i]` As for `(i)`, but without a bounds check. Not recommended.
- `(r,c)` For `Mat` and 2D field classes, access the element/object stored at row `r` and column `c`. An exception is thrown if the requested element is out of bounds.
- `.at(r,c)` As for `(r,c)`, but without a bounds check. Not recommended.
- `(r,c,s)` For `Cube` and 3D field classes, access the element/object stored at row `r`, column `c`, and slice `s`. An exception is thrown if the requested element is out of bounds.
- `.at(r,c,s)` As for `(r,c,s)`, but without a bounds check. Not recommended.

## 8.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> access1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++
      A(1,1) = 123.0;     // set element at row 2, column 2

      vec B(2, fill::randu);

      double x = A(0,1);  // copy element at row 1, column 2 to a double
      double y = B(1);    // copy element at coordinate 2 to a double

      uword i, j; // int also works
      uword N = A.n_rows;
      uword M = A.n_cols;

      for(i = 0; i < N; ++i) {
        for(j = 0; j < M; ++j) {
          A(i,j) = A(i,j) + x + y;
        }
      }

      return as_doubles_matrix(A);  // convert from C++ to R
    }

## 8.2 Caveats

For `.at()` or `[i]`, `.at(r,c)` and `.at(r,c,s)`:

- Indexing in C++ starts at 0
- Accessing elements without bounds checks is slightly faster, but is not recommended until your code has been thoroughly debugged first
- Accessing elements via `[r,c]` and `[r,c,s]` does not work correctly in C++; instead use `(r,c)` and `(r,c,s)`

The indices of elements are specified via the `uword` type, which is a `typedef` for an unsigned integer type. When using loops to access elements, it more efficient to use `uword` instead of `int`.

# 9 Element initialisation

Set elements in `Mat`, `Col` and `Row` via braced initialiser lists.

## 9.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> initialization1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++
      mat B = {{1, 2}, {3, 4}}; // create new matrix
      vec C = {1, 2}; // create new column vector

      // sum C to the diagonal of A
      A(0,0) = A(0,0) + C(0);
      A(1,1) = A(1,1) + C(1);

      mat D = A + B;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

# 10 Zeros

Set the elements of an object to zero, optionally first changing the size to specified dimensions.

`.zeros()` (member function of `Mat`, `Col`, `Row`, `SpMat`, `Cube`) `.zeros(n_elem)` (member function of `Col` and `Row`) `.zeros(n_rows, n_cols)` (member function of `Mat` and `SpMat`) `.zeros(n_rows, n_cols, n_slices)` (member function of `Cube`) `.zeros(size(X))` (member function of `Mat`, `Col`, `Row`, `Cube`, `SpMat`)

## 10.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> zeros1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++
      A.zeros();          // set all elements to zero

      mat B;
      B.zeros(size(A)); // set size to be the same as A and set all elements to zero

      mat C(A.n_rows, A.n_cols, fill::zeros);

      mat D = A + B + C;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

# 11 Ones

Set all the elements of an object to one, optionally first changing the size to specified dimensions.

Function | Mat | Col | Row | Cube
--- | --- | --- | --- | ---
`.ones()` | ✓ | ✓ | ✓ | ✓
`.ones(n_elem)` |  | ✓ | ✓ |
`.ones(n_rows, n_cols)` | ✓ |  |  |
`.ones(n_rows, n_cols, n_slices)` |  |  |  | ✓
`.ones(size(X))` | ✓ | ✓ | ✓ | ✓

## 11.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> ones1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++
      A.ones();          // set all elements to zero

      mat B;
      B.ones(size(A)); // set size to be the same as A and set all elements to zero

      mat C(A.n_rows, A.n_cols, fill::ones);

      mat D = A + B + C;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

# 12 Eye

`.eye()` is member function of `Mat` and `SpMat`. `.eye(n_rows, n_cols)` sets the elements along the main diagonal to one and off-diagonal elements to zero, optionally first changing the size to specified dimensions. `.eye(size(X))` creates an identity matrix is generated when `n_rows = n_cols`.

## 12.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> eye1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++
      A.eye();            // create an identity matrix

      mat B;
      B.eye(size(A)); // another identity matrix

      uword N = A.n_rows;
      uword M = A.n_cols;
      mat C(N, M, fill::randu);
      C.eye(N, M); // yet another identity matrix

      mat D = A + B + C;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

# 13 Random uniform

Set all the elements to random values from a uniform distribution in the [0,1] interval, optionally first changing the size to specified dimensions.

For complex elements, the real and imaginary parts are treated separately.

Function/Method | Mat | Col | Row | Cube
--- | --- | --- | --- | ---
`.randu()` | ✓ | ✓ | ✓ | ✓
`.randu(n_elem)` |  | ✓ | ✓ |
`.randu(n_rows, n_cols)` | ✓ |  |  |
`.randu(n_rows, n_cols, n_slices)` |  |  |  | ✓
`.randu(size(X))` | ✓ | ✓ | ✓ | ✓

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> randu1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++

      mat B;
      B.randu(size(A)); // random uniform matrix with the same size as A

      mat C(A.n_rows, A.n_cols, fill::randu);

      mat D = A + B + C;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> randu2_(const int& n) {
      GetRNGstate();  // Ensure R's RNG state is synchronized
      mat y(n, n);
      ::arma_rng::randu<double>::fill(y.memptr(), y.n_elem);
      PutRNGstate();

      return as_doubles_matrix(y);
    }

# 14 Normal distribution

Set all the elements to random values from a normal distribution with zero mean and unit variance, optionally first changing the size to specified dimensions.

For complex elements, the real and imaginary parts are treated separately.

Function/Method | Mat | Col | Row | Cube
--- | --- | --- | --- | ---
`.randn()` | ✓ | ✓ | ✓ | ✓
`.randn(n_elem)` |  | ✓ | ✓ |
`.randn(n_rows, n_cols)` | ✓ |  |  |
`.randn(n_rows, n_cols, n_slices)` |  |  |  | ✓
`.randn(size(X))` | ✓ | ✓ | ✓ | ✓

## 14.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> randn1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++

      mat B;
      B.randn(size(A)); // random normal matrix with the same size as A

      mat C(A.n_rows, A.n_cols, fill::randn);

      mat D = A + B + C;

      return as_doubles_matrix(D);  // convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> randn2_(const int& n) {
      GetRNGstate();  // Ensure R's RNG state is synchronized
      mat y(n, n);
      ::arma_rng::randn<double>::fill(y.memptr(), y.n_elem);
      PutRNGstate();

      return as_doubles_matrix(y);
    }

# 15 Fill

Sets the elements to a specified value

`.fill(value)` is a member function of `Mat`, `Col`, `Row`, `Cube`, `field`.

The type of value must match the type of elements used by the container object (e.g., for `Mat` the type is `double`)

## 15.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> fill1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++

      uword N = A.n_rows;
      uword M = A.n_cols;

      mat B(size(A), fill::value(200.0));  // create a matrix filled with 200.0
      mat C(N, M, fill::value(100.0));     // matrix filled with 100.0
      mat D(N, M, fill::zeros);            // matrix filled with zeros
      mat E(N, M, fill::ones);             // matrix filled with ones

      mat F = A + B + C + D + E;

      return as_doubles_matrix(F);  // convert from C++ to R
    }

# 16 Imbue

`.imbue(functor)` is a member function of `Mat`, `Col`, `Row` and `Cube`, it fills the elements with values provided by a functor. The argument can be a functor or lambda function.

For matrices, filling is done column-by-column (e.g., column 0 is filled, then column 1, etc.)

For cubes, filling is done slice-by-slice, with each slice treated as a matrix

## 16.1 Examples

    a <- matrix(runif(4), nrow = 2, ncol = 2)

    [[cpp4r::register]] doubles_matrix<> imbue1_(const doubles_matrix<>& a) {
      mat A = as_Mat(a);  // convert from R to C++

      std::mt19937 engine;  // Mersenne twister random number engine
      std::uniform_real_distribution<double> distr(0.0, 1.0);

      mat B(size(A), fill::none);                // create an empty matrix
      B.imbue([&]() { return distr(engine); });  // fill with random values

      mat C = A + B;

      return as_doubles_matrix(C);  // convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> imbue2_(const doubles_matrix<>& a) {
      GetRNGstate();  // Ensure R's RNG state is synchronized

      mat A = as_Mat(a);  // Convert from R to C++

      mat B(size(A), fill::none);  // Create an empty matrix
      B.imbue([]() { return unif_rand(); });  // Fill with random values

      mat C = A + B;

      PutRNGstate();

      return as_doubles_matrix(C);  // Convert from C++ to R
    }

# 17 Clean

`.clean(threshold)` is a member function of `Mat`, `Col`, `Row`, `Cube`, and `SpMat`. It can be used to sparsify a matrix, in the sense of zeroing values with small magnitudes.

- For objects with non-complex elements: each element with an absolute value less or equal to the threshold is replaced by zero.
- For objects with complex elements: for each element, each component (real and imaginary) with an absolute value less or equal to the threshold is replaced by zero.

## 17.1 Examples

    [[cpp4r::register]] doubles_matrix<> clean1_(const int& n) {
      mat A(n, n, fill::randu); // create a random matrix

      A(0, 0) = datum::eps; // set the diagonal with small values (+/- epsilon)
      A(1, 1) = -datum::eps;

      A.clean(datum::eps); // set elements with small values to zero

      return as_doubles_matrix(A); // Convert from C++ to R
    }

## 17.2 Caveat

To explicitly convert from dense storage to sparse storage, use the `SpMat`.

# 18 Replace

`.replace( old_value, new_value )` is a member function of `Mat`, `Col`, `Row`, `Cube`, and `SpMat`.

For all elements equal to `old_value`, set them to `new_value`.

## 18.1 Examples

    [[cpp4r::register]] doubles_matrix<> replace1_(const int& n) {
      mat A(n, n, fill::randu); // create a random matrix

      A.diag().fill(datum::nan); // set the diagonal with NaN values
      A.replace(datum::nan, 0);  // replace each NaN with 0

      return as_doubles_matrix(A); // Convert from C++ to R
    }

## 18.2 Caveats

- The type of `old_value` and `new_value` must match the type of elements used by the container object (e.g., for `Mat` the type is `double`).
- Floating point numbers (`float` and `double`) are approximations due to their limited precision.
- For sparse matrices (`SpMat`), replacement is not done when `old_value = 0`.

# 19 Clamp

`.clamp(min_value, max_value)` is a member function of `Mat`, `Col`, `Row`, `Cube` and `SpMat` that transforms all values lower than `min_val` to `min_val`, and all values higher than `max_val` to `max_val`.

- For complex elements, the real and imaginary components are clamped separately.
- For sparse matrices, clamping is applied only to the non-zero elements.

## 19.1 Examples

    [[cpp4r::register]] doubles_matrix<> clamp1_(const int& n) {
      mat A(n, n, fill::randu); // create a random matrix
      A.diag().fill(0.1);       // set the diagonal with 0.1 values

      A.clamp(0.2, 0.8); // clamp values to the [0.2, 0.8] interval

      return as_doubles_matrix(A); // Convert from C++ to R
    }

# 20 Transform

`.transform(functor)` is a member function of `Mat`, `Col`, `Row`, `Cube`, and `SpMat`. The argument can be a functor or lambda function.

- For dense matrices, transformation is done column-by-column for all elements.
- For sparse matrices, transformation is done column-by-column for non-zero elements.
- For cubes, transformation is done slice-by-slice, with each slice treated as a matrix.

## 20.1 Examples

    [[cpp4r::register]] doubles_matrix<> transform1_(const int& n) {
      mat A(n, n, fill::ones);  // create a matrix filled with ones
      A.transform([](double val) { return (val + 122.0); });
      return as_doubles_matrix(A); // Convert from C++ to R
    }

# 21 For each

`.for_each(functor)` is a member function of `Mat`, `Col`, `Row`, `Cube`, `SpMat`, and `field`. The argument can be a functor or lambda function.

- For dense matrices and fields, the processing is done column-by-column for all elements.
- For sparse matrices, the processing is done column-by-column for non-zero elements.
- For cubes, processing is done slice-by-slice, with each slice treated as a matrix.

## 21.1 Examples

    [[cpp4r::register]] doubles_matrix<> for_each1_(const int& n) {
      // add 122 to each element in a dense matrix, the '&' is important
      mat D(n, n, fill::ones);
      D.for_each([](mat::elem_type& val) { val += 122.0; });

      // add 122 to each non-zero element in a sparse matrix
      sp_mat S;
      S.sprandu(n, n, 1.0);
      S.for_each([](sp_mat::elem_type& val) { val += 123.0; });

      // set the size of all matrices in a field
      field<mat> F(2, 2);
      F.for_each([n](mat& X) { X.zeros(n, n); });  // capture n for the lambda

      mat res = D + S + F(0) + F(1);

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 22 Set size

Change the size of an object, without explicitly preserving data and without initialising the elements (e.g., elements may contain garbage values, including `NaN`).

- `.set_size(n_elem)` (member function of `Col`, `Row`, `field`)
- `.set_size(n_rows, n_cols)` (member function of `Mat`, `SpMat`, `field`)
- `.set_size(n_rows, n_cols, n_slices)` (member function of `Cube` and `field`)
- `.set_size(size(X))` (member function of `Mat`, `Col`, `Row`, `Cube`, `SpMat`, `field`)

To initialise the elements to zero while changing the size, use `.zeros()` instead. To explicitly preserve data while changing the size, use `.reshape()` or `.resize()` instead.

## 22.1 Examples

    [[cpp4r::register]] doubles set_size1_(const int& n) {
      mat A;
      A.set_size(n, n);  // or:  mat A(n, n, fill::none);

      mat B;
      B.set_size(size(A));  // or:  mat B(size(A), fill::none);

      vec C;
      C.set_size(n);  // or:  vec v(n, fill::none);

      A.fill(1.0);  // set all elements to 1.0
      B.fill(2.0);  // set all elements to 2.0
      C.fill(3.0);  // set all elements to 3.0

      vec res = A.col(0) + B.col(1) + C;

      return as_doubles(res);  // Convert from C++ to R
    }

# 23 Reshape

Recreate an object according to given size specifications, with the elements taken from the previous version of the object in a column-wise manner. The elements in the generated object are placed column-wise (e.g., the first column is filled up before filling the second column)

- `.reshape(n_rows, n_cols)` (member function of `Mat` and `SpMat`)
- `.reshape(n_rows, n_cols, n_slices)` (member function of `Cube`)
- `.reshape(size(X))` (member function of `Mat`, `Cube`, `SpMat`)

The layout of the elements in the recreated object will be different to the layout in the previous version of the object

If the total number of elements in the previous version of the object is less than the specified size, the extra elements in the recreated object are set to zero

If the total number of elements in the previous version of the object is greater than the specified size, only a subset of the elements is taken

## 23.1 Examples

    [[cpp4r::register]] doubles_matrix<> reshape1_(const int& n) {
      mat A(n + 1, n - 1, fill::randu);
      A.reshape(n - 1, n + 1);
      return as_doubles_matrix(A);  // Convert from C++ to R
    }

## 23.2 Caveats

- `.reshape()` is considerably slower than `.set_size()`.
- To change the size without preserving data, use `.set_size()`.
- To grow/shrink the object while preserving the elements and the layout of the elements, use `.resize()`
- To flatten a matrix into a vector, use `vectorise()` or `.as_col()`/`.as_row()`.

# 24 Resize

Resize an object according to given size specifications, while preserving the elements and the layout of the elements. It can be used for growing or shrinking an object (e.g., adding/removing rows, and/or columns, and/or slices).

- `.resize(n_elem)`: member function of `Col`, `Row`.
- `.resize(n_rows, n_cols)`: member function of `Mat` and `SpMat`.
- `.resize(n_rows, n_cols, n_slices)`: member function of `Cube`.
- `.resize(size(X))`: member function of `Mat`, `Col`, `Row`, `Cube`, `SpMat`.

## 24.1 Examples

    [[cpp4r::register]] doubles_matrix<> resize1_(const int& n) {
      mat A(n + 1, n - 1, fill::randu);
      A.resize(n - 1, n + 1);
      return as_doubles_matrix(A);  // Convert from C++ to R
    }

## 24.2 Caveats

- `.resize()` is considerably slower than `.set_size()`.
- to change the size without preserving data, `.set_size()` instead.

# 25 Copy size

`.copy_size(A)` sets the size of a matrix/vector/cube to be the same as matrix/vector/cube `A`.

## 25.1 Examples

    [[cpp4r::register]] integers copy_size1_(const int& n) {
      mat A(n, n, fill::randu);

      mat B;
      B.copy_size(A);

      int N = B.n_rows;
      int M = B.n_cols;

      writable::integers res({N, M});
      res.attr("names") = strings({"n_rows", "n_cols"});

      return as_integers(res);  // Convert from C++ to R
    }

## 25.2 Caveat

To set the size of an object `B`, `A` must be of the same type as `B`. For example, the size of a matrix cannot be set by providing a cube.

# 26 Reset

`.reset()` sets a matrix/vector size to zero (the object will have no elements).

## 26.1 Examples

    [[cpp4r::register]] integers reset1_(const int& n) {
      mat A(n, n, fill::randu);
      A.reset();

      int N = A.n_rows;
      int M = A.n_cols;

      writable::integers res({N, M});
      res.attr("names") = strings({"n_rows", "n_cols"});

      return as_integers(res);  // Convert from C++ to R
    }

# 27 Submatrix views

A collection of member functions of `Mat`, `Col` and `Row` classes that provide read/write access to submatrix views.

## 27.1 Contiguous views for matrix

- `X.col(col_number)`
- `X.row(row_number)`
- `X.cols(first_col, last_col)`
- `X.rows(first_row, last_row)`
- `X.submat(first_row, first_col, last_row, last_col)`
- `X(span(first_row, last_row), span(first_col, last_col))`
- `X(first_row, first_col, size(n_rows, n_cols))`
- `X(first_row, first_col, size(Y))` (`Y` is a matrix)
- `X(span(first_row, last_row), col_number)`
- `X(row_number, span(first_col, last_col))`
- `X.head_cols(number_of_cols)`
- `X.head_rows(number_of_rows)`
- `X.tail_cols(number_of_cols)`
- `X.tail_rows(number_of_rows)`
- `X.unsafe_col(col_number)` (use with caution)

## 27.2 Contiguous views for vector

- `Y(span(first_index, last_index))`
- `Y.subvec(first_index, last_index)`
- `Y.subvec(first_index, size(X))` (`X` is a vector)
- `Y.head(number_of_elements)`
- `Y.tail(number_of_elements)`

## 27.3 Non-contiguous views for matrix or vector:

- `X.elem(vector_of_indices)`
- `X(vector_of_indices)`
- `X.cols(vector_of_column_indices)`
- `X.rows(vector_of_row_indices)`
- `X.submat(vector_of_row_indices, vector_of_column_indices)`
- `X(vector_of_row_indices, vector_of_column_indices)`

Instances of `span(start, end)` can be replaced by `span::all_` to indicate the entire range.

For functions requiring one or more vector of indices, for example `X.submat(vector_of_row_indices, vector_of_column_indices)`, each vector of indices must be of type `uvec`.

In the function `X.elem(vector_of_indices)`, elements specified in `vector_of_indices` are accessed. `X` is interpreted as one long vector, with column-by-column ordering of the elements of `X`. The `vector_of_indices` must evaluate to a vector of type `uvec` (e.g., generated by the `find()` function). The aggregate set of the specified elements is treated as a column vector (e.g., the output of `X.elem()` is always a column vector).

The function `.unsafe_col()` is provided for speed reasons and should be used only if you know what you are doing. It creates a seemingly independent `Col` vector object (e.g., `vec`), but uses memory from the existing matrix object. As such, the created vector is not alias safe, and does not take into account that the underlying matrix memory could be freed (e.g., due to any operation involving a size change of the matrix).

## 27.4 Examples

    [[cpp4r::register]] doubles_matrix<> subview1_(const int& n) {
      mat A(n, n, fill::zeros);

      A.submat(0,1,2,3) = randu<mat>(3,3);
      A(span(0,2), span(1,3)) = randu<mat>(3,3);
      A(0,1, size(3,3)) = randu<mat>(3,3);

      mat B = A.submat(0,1,2,3);
      mat C = A(span(0,2), span(1,3) );
      mat D = A(0, 1, size(3,3) );

      A.col(1) = randu<mat>(5,1);
      A(span::all, 1) = randu<mat>(5,1);

      mat X(5, 5, fill::randu);

      // get all elements of X that are greater than 0.5
      vec q = X.elem( find(X > 0.5) );

      // add 123 to all elements of X greater than 0.5
      X.elem( find(X > 0.5) ) += 123.0;

      // set four specific elements of X to 1
      uvec indices = { 2, 3, 6, 8 };

      X.elem(indices) = ones<vec>(4);

      // add 123 to the last 5 elements of vector a
      vec a(10, fill::randu);
      a.tail(5) += 123.0;

      // add 123 to the first 3 elements of column 2 of X
      X.col(2).head(3) += 123;

      return as_doubles_matrix(X);  // Convert from C++ to R
    }

# 28 Subcube views and slices

A collection of member functions of the `Cube` class that provide subcube views.

## 28.1 Contiguous views for cube

- `Q.slice(slice_number)`
- `Q.slices(first_slice, last_slice)`
- `Q.row(row_number)`
- `Q.rows(first_row, last_row)`
- `Q.col(col_number)`
- `Q.cols(first_col, last_col)`
- `Q.subcube( first_row, first_col, first_slice, last_row, last_col, last_slice)`
- `Q(span(first_row, last_row), span(first_col, last_col), span(first_slice, last_slice))`
- `Q(first_row, first_col, first_slice, size(n_rows, n_cols, n_slices))`
- `Q(first_row, first_col, first_slice, size(R))` (`R` is a cube)
- `Q.head_slices(number_of_slices)`
- `Q.tail_slices(number_of_slices)`
- `Q.tube(row, col)`
- `Q.tube(first_row, first_col, last_row, last_col)`
- `Q.tube(span(first_row, last_row), span(first_col, last_col))`
- `Q.tube(first_row, first_col, size(n_rows, n_cols))`

## 28.2 Non-contiguous views for cube

`Q.elem(vector_of_indices)`, `Q(vector_of_indices)`, and `Q.slices( vector_of_slice_indices)` are instances of `span(a,b)` that can be replaced by:

- `span()` or `span::all`, to indicate the entire range.
- `span(a)`, to indicate a particular row, column or slice.

An individual slice, accessed via `.slice()`, is an instance of the `Mat` class (a reference to a matrix is provided).

All `.tube()` forms are variants of `.subcube()`, using `first_slice = 0` and `last_slice = Q.n_slices-1`. The `.tube(row,col)` form uses `row = first_row = last_row`, and `col = first_col = last_col`.

In the function `Q.elem(vector_of_indices)`, elements specified in `vector_of_indices` are accessed. `Q` is interpreted as one long vector, with slice-by-slice and column-by-column ordering of the elements of `Q`. The `vector_of_indices` must evaluate to a vector of type `uvec` (e.g., generated by the `find()` function). The aggregate set of the specified elements is treated as a column vector (e.g., the output of `Q.elem()` is always a column vector).

In the function `Q.slices(vector_of_slice_indices)`, slices specified in `vector_of_slice_indices` are accessed. The `vector_of_slice_indices` must evaluate to a vector of type `uvec`.

## 28.3 Examples

    [[cpp4r::register]] doubles_matrix<> subview2_(const int& n) {
      cube A(n, 3, 4, fill::randu);

      mat B = A.slice(1); // each slice is a matrix

      A.slice(0) = randu<mat>(2,3);
      A.slice(0)(1,2) = 99.0;

      A.subcube(0,0,1,  1,1,2)           = randu<cube>(2,2,2);
      A(span(0,1), span(0,1), span(1,2)) = randu<cube>(2,2,2);
      A(0,0,1, size(2,2,2))              = randu<cube>(2,2,2);

      // add 123 to all elements of A greater than 0.5
      A.elem( find(A > 0.5) ) += 123.0;

      cube C = A.head_slices(2);  // get first two slices

      A.head_slices(2) += 123.0;

      mat res = A.slice(0) + B + C.slice(1);

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 29 Subfield views

A collection of member functions of the `field` class that provide subfield views.

For a 2D field `F`, the subfields are accessed as:

- `F.row(row_number)`
- `F.col(col_number)`
- `F.rows(first_row, last_row)`
- `F.cols(first_col, last_col)`
- `F.subfield(first_row, first_col, last_row, last_col)`
- `F(span(first_row, last_row), span(first_col, last_col))`
- `F(first_row, first_col, size(G))` (`G` is a 2D field)
- `F(first_row, first_col, size(n_rows, n_cols))`

For a 3D field `F`, the subfields are accessed as:

- `F.slice(slice_number)`
- `F.slices(first_slice, last_slice)`
- `F.subfield(first_row, first_col, first_slice, last_row, last_col, last_slice)`
- `F(span(first_row, last_row), span(first_col, last_col), span(first_slice, last_slice))`
- `F(first_row, first_col, first_slice, size(G))` (`G` is a 3D field)
- `F(first_row, first_col, first_slice, size(n_rows, n_cols, n_slices))`

Instances of `span(a,b)` can be replaced by:

- `span()` or `span::all`, to indicate the entire range.
- `span(a)`, to indicate a particular row or column.

# 30 Diagonal

`.diag()` is a member functions of `Mat` and `SpMat` with read/write access to the diagonal in a matrix. The argument can be empty or a value `k` to specify the diagonal to (`k = 0` by default). The diagonal is interpreted as a column vector within expressions.

- `k = 0` indicates the main diagonal (default setting)
- `k < 0` indicates the `k`-th sub-diagonal (below main diagonal, towards bottom-left corner)
- `k > 0` indicates the `k`-th super-diagonal (above main diagonal, towards top-right corner)

## 30.1 Examples

    [[cpp4r::register]] doubles diagonal1_(const int& n) {
      mat X(n, n, fill::randu);

      vec A = X.diag(); // extract the main diagonal
      double B = accu(X.diag(1)); // sum of elements on the first upper diagonal
      double C = accu(X.diag(-1)); // sum of elements on the first lower diagonal

      X.diag() = randu<vec>(n);
      X.diag() += A;
      X.diag() /= B;
      X.diag() *= C;

      sp_mat S = sprandu<sp_mat>(n, n, 0.0);
      S.diag().ones();

      vec v(S.diag());  // copy sparse diagonal to dense vector
      v += X.diag();

      return as_doubles(v);  // Convert from C++ to R
    }

## 30.2 Caveat

To calculate only the diagonal elements of a compound expression, use `diagvec()` or `diagmat()`.

# 31 Each col

`.each_col()` is a member function of `Mat`. It applies a vector operation to each column of a matrix, and are similar to “broadcasting” in Matlab/Octave. The argument can be empty, a vector of indices, or a lambda function.

Operation | `.each_col()` | `.each_col(vector_of_indices)` | `.each_col(lambda)`
--- | --- | --- | ---
`+` addition | ✓ | ✓ |
`+=` in-place addition | ✓ | ✓ |
`-` subtraction | ✓ | ✓ |
`-=` in-place subtraction | ✓ | ✓ |
`%` element-wise multiplication | ✓ | ✓ |
`%=` in-place element-wise multiplication | ✓ | ✓ |
`/` element-wise division | ✓ | ✓ |
`/=` in-place element-wise division | ✓ | ✓ |
`=` assignment (copy) | ✓ | ✓ |
`lambda` (lambda function) |  |  | ✓

## 31.1 Examples

    [[cpp4r::register]] doubles_matrix<> each_col1_(const int& n) {
      mat X(n, n + 1, fill::ones);

      // create a vector with n elements ranging from 5 to 10
      vec v = linspace<vec>(5, 10, n);

      // in-place addition of v to each column vector of X
      X.each_col() += v;

      // generate Y by adding v to each column vector of X
      mat Y = X.each_col() + v;

      // subtract v from columns 1 and 2 of X
      X.cols(0, 1).each_col() -= v;

      uvec indices(2);
      indices(0) = 1;
      indices(1) = 2;

      X.each_col(indices) = v;  // copy v to columns 1 and 2 of X

      // lambda function with non-const vector
      X.each_col([](vec& a) { 2 * a; });

      const mat& XX = X;

      // lambda function with const vector
      XX.each_col([](const vec& b) { 3 * b; });

      mat res = X + Y + XX;

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 32 Each row

`.each_row()`, `.each_row(vector_of_indices)`, `.each_row(lambdaction)` are member functions of `Mat`. These apply a vector operation to each row of a matrix, and are similar to “broadcasting” in Matlab/Octave.

## 32.1 Form 1

`.each_row()` supports the following operations:

- `+` addition
- `+=` in-place addition
- `-` subtraction
- `-=` in-place subtraction
- `%` element-wise multiplication
- `%=` in-place element-wise multiplication
- `/` element-wise division
- `/=` in-place element-wise division
- `=` assignment (copy)

## 32.2 Form 2

`.each_row(vector_of_indices)` supports the same operations as form 1. The argument `vector_of_indices` contains a list of indices of the rows to be used, and it must evaluate to a vector of type `uvec`.

## 32.3 Form 3

`.each_col(lambdaction)` applies the given `lambdaction` to each column vector. The function must accept a reference to a `Row` object with the same element type as the underlying matrix.

## 32.4 Examples

    [[cpp4r::register]] doubles_matrix<> each_row1_(const int& n) {
      mat X(n + 1, n, fill::ones);

      // create a vector with n elements ranging from 5 to 10
      rowvec v = linspace<rowvec>(5, 10, n);

      // in-place addition of v to each rows vector of X
      X.each_row() += v;

      // generate Y by adding v to each rows vector of X
      mat Y = X.each_row() + v;

      // subtract v from rows 1 and 2 of X
      X.rows(0, 1).each_row() -= v;

      uvec indices(2);
      indices(0) = 1;
      indices(1) = 2;

      X.each_row(indices) = v;       // copy v to columns 1 and 2 of X

      // lambda function with non-const vector
      X.each_row([](rowvec& a) { a / 2; });

      const mat& XX = X;

      // lambda function with const vector
      XX.each_row([](const rowvec& b) { b / 3; });

      mat res = X + Y + XX;

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 33 Each slice

`.each_slice()` is a member function of `Cube` that applies a matrix operation to each slice of a cube, with each slice treated as a matrix. It is similar to “broadcasting” in Matlab/Octave.

## 33.1 Form 1

`.each_slice(vector_of_indices)`

Supported operations:

- `+` addition
- `+=` in-place addition
- `-` subtraction
- `-=` in-place subtraction
- `%` element-wise multiplication
- `%=` in-place element-wise multiplication
- `/` element-wise division
- `/=` in-place element-wise division
- `*` matrix multiplication
- `*=` in-place matrix multiplication
- `=` assignment (copy)

## 33.2 Form 2

`.each_slice(lambdaction)`

- The argument vector_of_indices contains a list of indices of the slices to be used; it must evaluate to a vector of type `uvec`.
- Arithmetic operations as per form 1 are supported, except for `*` and `*=` (e.g., matrix multiplication).

## 33.3 Form 3

`.each_slice(lambdaction, use_mp)`

- Apply the given `lambdaction` to each slice.
- The function must accept a reference to a `Mat` object with the same element type as the underlying cube.

## 33.4 Form 4

- Apply the given `lambdaction` to each slice, as per form 3.
- The argument `use_mp` is a bool to enable the use of OpenMP for multi-threaded execution of `lambdaction` on multiple slices at the same time.
- The order of processing the slices is not deterministic (e.g., slice 2 can be processed before slice 1).
- `lambdaction` must be thread-safe, e.g., it must not write to variables outside of its scope.

## 33.5 Examples:

    [[cpp4r::register]] doubles_matrix<> each_slice1_(const int& n) {
      cube C(n, n + 1, 6, fill::randu);

      mat M = repmat(linspace<vec>(1, n, n), 1, n + 1);

      C.each_slice() += M;  // in-place addition of M to each slice of C

      cube D = C.each_slice() + M;  // generate D by adding M to each slice of C

      // sum all slices of D into a single n x (n + 1) matrix
      mat D_flat = sum(D, 2);

      uvec indices(2);
      indices(0) = 2;
      indices(1) = 4;

      C.each_slice(indices) = M;  // copy M to slices 2 and 4 in C
      C.each_slice([](mat& X) { X * 2.0; });  // lambda function with non-const matrix
      mat C_flat = sum(C, 2);

      const cube& CC = C;
      CC.each_slice([](const mat& X) { X / 3.0; });  // lambda function with const matrix

      mat CC_flat = sum(CC, 2);

      mat res = C_flat + D_flat + CC_flat;

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 34 Set real

`.set_real(X)` sets the real part of an object. `X` must have the same size as the recipient object.

## 34.1 Examples

    [[cpp4r::register]] complexes_matrix<> set_real1_(const int& n) {
      mat A(n + 1, n - 1, fill::randu);

      cx_mat C(n + 1, n - 1, fill::zeros);

      C.set_real(A);

      return as_complexes_matrix(C);  // Convert from C++ to R
    }

## 34.2 Caveat

To directly construct a complex matrix out of two real matrices, the following code is faster:

    [[cpp4r::register]] complexes_matrix<> set_real2_(const int& n) {
      mat A(n - 1, n + 1, fill::randu);
      mat B(n - 1, n + 1, fill::randu);

      cx_mat C = cx_mat(A,B);

      return as_complexes_matrix(C);  // Convert from C++ to R
    }

# 35 Set imaginary

`.set_imaginary(X)` sets the imaginary part of an object. `X` must have the same size as the recipient object.

## 35.1 Examples

    [[cpp4r::register]] complexes_matrix<> set_imag1_(const int& n) {
      mat B(n + 1, n - 1, fill::randu);

      cx_mat C(n + 1, n - 1, fill::zeros);

      C.set_imag(B);

      return as_complexes_matrix(C);  // Convert from C++ to R
    }

## 35.2 Caveat

To directly construct a complex matrix out of two real matrices, the following code is faster:

    [[cpp4r::register]] complexes_matrix<> set_imag2_(const int& n) {
      mat A(n - 1, n + 1, fill::randu);
      mat B(n - 1, n + 1, fill::randu);

      cx_mat C = cx_mat(A,B);

      return as_complexes_matrix(C);  // Convert from C++ to R
    }

# 36 Insert columns

`.insert_cols()` is a member function of `Mat`, `Row` and `Cube`. The arguments can be `colnumber, X` to indicate the column number and the matrix to insert, or `colnumber, number_of_cols` to indicate the column number and the number of columns to insert.

The `X` argument inserts a copy of `X` at the specified column. `X` must have the same number of rows (and slices) as the recipient object.

The `number_of_cols` argument expands the object by creating new columns that are set to zero.

## 36.1 Examples

    [[cpp4r::register]] doubles_matrix<> insert_columns1_(const int& n) {
      mat A(n, n * 2, fill::randu);
      mat B(n, n - 1, fill::ones);

      // at column n - 1, insert a copy of B
      // A will now have 3n - 1 columns
      A.insert_cols(n - 1, B);

      // at column 1, insert 2n zeroed columns
      // B will now have 3n - 1 columns
      B.insert_cols(1, n * 2);

      mat res = A + B;

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 37 Insert rows

`.insert_rows()` is a member function of `Mat`, `Row` and `Cube`. The arguments can be `rownumber, X` to indicate the row number and the matrix to insert, or `rownumber, number_of_rows` to indicate the row number and the number of rows to insert.

The `X` argument inserts a copy of `X` at the specified column. `X` must have the same number of columns (and slices) as the recipient object.

The `number_of_rows` argument expands the object by creating new rows that are set to zero.

## 37.1 Examples

    [[cpp4r::register]] doubles_matrix<> insert_rows1_(const int& n) {
      mat A(n * 2, n, fill::randu);
      mat B(n - 1, n, fill::ones);

      // at row n - 1, insert a copy of B
      // A will now have 3n - 1 rows
      A.insert_rows(n - 1, B);

      // at row 1, insert 2n zeroed rows
      // B will now have 3n - 1 columns
      B.insert_rows(1, n * 2);

      mat res = A + B;

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 38 Insert slice

`.insert_slices()` is a member function of `Cube`. The arguments can be `slice_number, X` to indicate the slice number and the matrix to insert, or `slice_number, number_of_slices` to indicate the slice number and the number of slices to insert.

The `X` argument inserts a copy of `X` at the specified slice. `X` must have the same number of columns and rows as the recipient object.

The `number_of_slices` argument expands the object by creating new slices that are set to zero.

## 38.1 Examples

    [[cpp4r::register]] doubles_matrix<> insert_slices1_(const int& n) {
      cube A(n, n, n * 2, fill::randu);
      cube B(n, n, n - 1, fill::ones);

      // At slice n - 1, insert a copy of B
      // A will now have 3n - 1 slices
      A.insert_slices(n - 1, B);

      // At slice 1, insert 2n zeroed slices
      // B will now have 3n - 1 slices
      B.insert_slices(1, n * 2);

      mat res = sum(A + B);

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 39 Shed columns

`.shed_col(row_number)` and `.shed_cols(first_row, last_row)` are member functions of `Mat`, `Col`, `SpMat`, and `Cube`. With a single scalar argument it remove the specified column, and with two scalar arguments it removes the specified range of columns.

`.shed_cols(vector_of_indices)` is a member function of `Mat` and `Col`. With a vector of indices it must evaluate to a vector of type `uvec` containing the indices of the columns to remove.

## 39.1 Examples

    [[cpp4r::register]] doubles_matrix<> shed_columns1_(const int& n) {
      mat A(n, n * 5, fill::randu);

      // remove the first column
      A.shed_col(0);

      // remove columns 1 and 2
      A.shed_cols(0, 1);

      // remove columns 2 and 4
      uvec indices(2);
      indices(0) = 1;
      indices(1) = 3;
      A.shed_cols(indices);

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

# 40 Shed rows

`.shed_row(row_number)` and `.shed_rows(first_row, last_row)` are member functions of `Mat`, `Col`, `SpMat`, and `Cube`. With a single scalar argument it remove the specified rows, and with two scalar arguments it removes the specified range of rows.

`.shed_rows(vector_of_indices)` is a member function of `Mat` and `Row`. With a vector of indices it must evaluate to a vector of type `uvec` containing the indices of the rows to remove.

## 40.1 Examples

    [[cpp4r::register]] doubles_matrix<> shed_rows1_(const int& n) {
      mat A(n * 5, n, fill::randu);

      // remove the first row
      A.shed_row(0);

      // remove rows 1 and 2
      A.shed_rows(0, 1);

      // remove rows 2 and 4
      uvec indices(2);
      indices(0) = 1;
      indices(1) = 3;
      A.shed_rows(indices);

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

# 41 Shed slices

`.shed_slices()` is a member function of `Cube`. With a single scalar argument it remove the specified slices, and with two scalar arguments it removes the specified range of slices. With a vector of indices it must evaluate to a vector of type `uvec` containing the indices of the rows to remove. The arguments can be `slice_number` to indicate the slice number to remove, `first_slice, last_slice` to indicate the range of slices to remove, or `vector_of_indices` to indicate the indices of the slices to remove.

## 41.1 Examples

    [[cpp4r::register]] doubles_matrix<> shed_slices1_(const int& n) {
      cube A(n, n, n * 5, fill::randu);

      // remove the first slice
      A.shed_slice(0);

      // remove slices 1 and 2
      A.shed_slices(0, 1);

      // remove slices 2 and 4
      uvec indices(2);
      indices(0) = 1;
      indices(1) = 3;
      A.shed_slices(indices);

      mat res = sum(A, 2);

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

# 42 Swap columns

`.swap_cols( col1, col2 )` is a member functions of `Mat`, `Col`, `Row`, and `SpMat`. It swaps the contents of the specified columns.

## 42.1 Examples

    [[cpp4r::register]] doubles_matrix<> swap_columns1_(const int& n) {
      mat A(n, n * 5, fill::randu);

      // swap columns 1 and 2
      A.swap_cols(0, 1);

      // swap columns 2 and 4
      A.swap_cols(1, 3);

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

# 43 Swap rows

`.swap_rows( col1, col2 )` is a member functions of `Mat`, `Col`, `Row`, and `SpMat`. It swaps the contents of the specified rows.

## 43.1 Examples

    [[cpp4r::register]] doubles_matrix<> swap_rows1_(const int& n) {
      mat A(n * 5, n, fill::randu);

      // swap rows 1 and 2
      A.swap_rows(0, 1);

      // swap rows 2 and 4
      A.swap_rows(1, 3);

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

# 44 Swap

`.swap( X )` is a member function of `Mat`, `Col`, `Row`, and `Cube`. It swaps the contents with object `X`.

## 44.1 Examples

    [[cpp4r::register]] doubles_matrix<> swap1_(const int& n) {
      mat A(n, n + 1, fill::zeros);
      mat B(n * 2, n - 1, fill::ones);

      A.swap(B);

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

# 45 Memory pointer

`.memptr()` is a member function of `Mat`, `Col`, `Row`, and `Cube`. It obtains a raw pointer to the memory used for storing elements. Data for matrices is stored in a column-by-column order. Data for cubes is stored in a slice-by-slice (matrix-by-matrix) order.

## 45.1 Examples

    [[cpp4r::register]] doubles_matrix<> memptr1_(const int& n) {
      mat A(n, n, fill::randu);
      const mat B(n, n, fill::randu);

      double* A_mem = A.memptr();
      const double* B_mem = B.memptr();

      // alter A_mem
      // B_mem is const, so it cannot be altered
      for (int i = 0; i < n * n; ++i) {
        A_mem[i] += 123.0 + B_mem[i];
      }

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

## 45.2 Caveats

- The pointer becomes invalid after any operation involving a size change or aliasing.
- This function is not recommended for use unless you know what you are doing.

# 46 Column pointer

`.colptr( col_number )` is a member function of the `Mat` class that obtains a raw pointer to the memory used by elements in the specified column.

## 46.1 Examples

    [[cpp4r::register]] doubles_matrix<> colptr1_(const int& n) {
      mat A(n, n, fill::randu);

      // pointer to the memory of the first column of A
      double* Acol1_mem = A.colptr(0);

      // alter memory
      for (int i = 0; i < n; ++i) {
        Acol1_mem[i] += 123.0;
      }

      return as_doubles_matrix(A);  // Convert from C++ to R
    }

## 46.2 Caveats

- The pointer becomes invalid after any operation involving a size change or aliasing.
- This function is not recommended for use unless you know what you are doing.
- It is safer to use `submat()` instead.

# 47 Iterators

Iterators for traverse over all elements within the specified range. These return the column/row/slice of an object as a `uword` type.

## 47.1 Member functions

Dense matrices and vectors (`Mat`, `Col`, and `Row`):

- `.begin()` is an iterator referring to the first element.
- `.end()` is an iterator referring to the past the end element.
- `.begin_col(col_number)` is an iterator referring to the first element of the specified column.
- `.end_col(col_number)` is an iterator referring to the past-the-end element of the specified column.
- `begin_row(row_number)` is an iterator referring to the first element of the specified row.
- `end_row(row_number)` is an iterator referring to the past-the-end element of the specified row.

Cubes (`Cube`):

- `begin()` is an iterator referring to the first element.
- `end()` is an iterator referring to the past-the-end element.
- `begin_slice(slice_number)` iterator referring to the first element of the specified slice.
- `end_slice(slice_number)` iterator referring to the past-the-end element of the specified slice.

Sparse matrices (`SpMat`):

- `begin()` is an iterator referring to the first element.
- `end()` is an iterator referring to the past-the-end element.
- `begin_col(col_number)` is an iterator referring to the first element of the specified column.
- `end_col(col_number)` is an iterator referring to the past-the-end element of the specified column.
- `begin_row(row_number)` is an iterator referring to the first element of the specified row.
- `end_row(row_number)` is an iterator referring to the past-the-end element of the specified row.

Dense submatrices and subcubes (`submatrix` and `subcube`):

- `span(row, col)` and `span(row, col, slice)` can be used to specify the range of elements to iterate over.

## 47.2 Iterator types

Dense matrices and vectors (`Mat`, `Col`, and `Row`):

- `mat::iterator`, `vec::iterator` and `rowvec::iterator` are random access iterators, for read/write access to elements (which are stored column by column).
- `mat::const_iterator`, `vec::const_iterator` and `rowvec::const_iterator` are random access iterators, for read-only access to elements (which are stored column by column)
- `mat::col_iterator`, `vec::col_iterator` and `rowvec::col_iterator` random access iterators, for read/write access to the elements of specified columns.
- `mat::const_col_iterator`, `vec::const_col_iterator` and `rowvec::const_col_iterator` are random access iterators, for read-only access to the elements of specified columns.
- `mat::row_iterator` is a bidirectional iterator, for read/write access to the elements of specified rows.
- `mat::const_row_iterator` is a bidirectional iterator, for read-only access to the elements of specified rows.
- `vec::row_iterator` and `rowvec::row_iterator` are random access iterators, for read/write access to the elements of specified rows.
- `vec::const_row_iterator` and `rowvec::const_row_iterator` are random access iterators, for read-only access to the elements of specified rows.

Cubes (`Cube`):

- `cube::iterator` is a random access iterator, for read/write access to elements. The elements are ordered slice by slice; the elements within each slice are ordered column by column.
- `cube::const_iterator` is a random access iterator, for read-only access to elements.
- `cube::slice_iterator` is a random access iterator, for read/write access to the elements of a particular slice. The elements are ordered column by column.
- `cube::const_slice_iterator` is a random access iterator, for read-only access to the elements of a particular slice.

Sparse matrices (`SpMat`):

- `sp_mat::iterator` is a bidirectional iterator, for read/write access to elements (which are stored column by column).
- `sp_mat::const_iterator` is a bidirectional iterator, for read-only access to elements (which are stored column by column).
- `sp_mat::col_iterator` is a bidirectional iterator, for read/write access to the elements of a specific column.
- `sp_mat::const_col_iterator` is a bidirectional iterator, for read-only access to the elements of a specific column.
- `sp_mat::row_iterator` is a bidirectional iterator, for read/write access to the elements of a specific row.
- `sp_mat::const_row_iterator` is a bidirectional iterator, for read-only access to the elements of a specific row.

## 47.3 Examples

    [[cpp4r::register]] doubles_matrix<> iterators1_(const int& n) {
      mat X(n, n + 1, fill::randu);

      mat::iterator it = X.begin();
      mat::iterator it_end = X.end();

      for (; it != it_end; ++it) {
        (*it) += 123.0;
      }

      mat::col_iterator col_it = X.begin_col(1);    // start of column 1
      mat::col_iterator col_it_end = X.end_col(n);  //   end of column n

      for (; col_it != col_it_end; ++col_it) {
        (*col_it) = 321.0;
      }

      return as_doubles_matrix(X);  // Convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> iterators2_(const int& n) {
      cube X(n, n + 1, n + 2, fill::randu);

      cube::iterator it = X.begin();
      cube::iterator it_end = X.end();

      for (; it != it_end; ++it) {
        (*it) += 123.0;
      }

      cube::slice_iterator s_it = X.begin_slice(1);    // start of slice 1
      cube::slice_iterator s_it_end = X.end_slice(n);  // end of slice n

      for (; s_it != s_it_end; ++s_it) {
        (*s_it) = 321.0;
      }

      mat res = sum(X, 2);

      return as_doubles_matrix(res);  // Convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> iterators3_(const int& n) {
      sp_mat X = sprandu<sp_mat>(n, n * 2, 0.1);

      sp_mat::iterator it = X.begin();
      sp_mat::iterator it_end = X.end();

      for (; it != it_end; ++it) {
        (*it) += 123.0;
      }

      return as_doubles_matrix(X);  // Convert from C++ to R
    }

    [[cpp4r::register]] doubles_matrix<> iterators4_(const int& n) {
      mat X(n, n, fill::randu);

      for (double& val : X(span(0, 1), span(1, 1))) {
        val = 123.0;
      }

      return as_doubles_matrix(X);  // Convert from C++ to R
    }

## 47.4 Caveats

- Writing a zero value into a sparse matrix through an iterator will invalidate all current iterators associated with the sparse matrix.
- To modify the non-zero elements in a safer manner, use `.transform()` or `.for_each()` instead of iterators.
- For `submatrix` and `subcube` the iterators are intended only to be used with range-based for loops. Any other use is not supported. For example, the direct use of the `.begin()` and `.end()` functions, as well as the underlying iterators types is not supported. The implementation of submatrices and subcubes uses short-lived temporary objects that are subject to automatic deletion, and as such are error-prone to handle manually.

# 48 Compatibility container functions

Member functions for the `Col` and `Row` classes to mimic the functionality of containers in the C++ standard library:

- `.front()` accesses the first element in a vector.
- `.back()` accesses the last element in a vector.

Member functions for the `Col`, `Row`, `Mat`, `Cube` and `SpMat` classes to mimic the functionality of containers in the C++ standard library:

- `.clear()` removes the elements from an object.
- `.empty()` returns `true` if the object has no elements and `false` if the object has one or more elements.
- `.size()` returns the total number of elements in an object.

## 48.1 Examples

    [[cpp4r::register]] doubles compatibility1_(const int& n) {
      vec X(n, fill::randu);

      writable::doubles res = {X.front(), X.back()};

      res.attr("names") = strings({"front", "back"});

      return res;
    }

    [[cpp4r::register]] integers compatibility2_(const int& n) {
      mat X(n, n, fill::randu);

      writable::integers res(2);
      res[0] = X.n_rows;

      X.clear();
      res[1] = X.n_rows;

      res.attr("names") = strings({"before", "after"});

      return res;
    }

# 49 Convert matrix to column

`.as_col()` is a member function of the `Mat` class, it returns a flattened version of the matrix as a column vector. Flattening is done by concatenating all columns.

## 49.1 Examples

    [[cpp4r::register]] doubles as_col1_(const int& n) {
      mat M(n, n + 1, fill::randu);
      vec V = M.as_col();
      return as_doubles(V);
    }

# 50 Convert matrix to row

`.as_row()` is a member function of the `Mat` class, it returns a flattened version of the matrix as a row vector. Flattening is done by concatenating all rows.

## 50.1 Examples

    [[cpp4r::register]] doubles as_row1_(const int& n) {
      mat M(n, n + 1, fill::randu);
      vec V = M.as_row();
      return as_doubles(V);
    }

## 50.2 Caveat

Converting columns to rows is faster than converting rows to columns.

# 51 Convert column to matrix

`.col_as_mat(col_number)` is a member function of the `Cube` class, it returns a matrix of the specified cube column and the number of rows is preserved. Given a cube of size `R x C x S`, the resultant matrix size is `R x S`.

## 51.1 Examples

    [[cpp4r::register]] list col_as_mat1_(const int& n) {
      cube C(n, n + 1, n + 2, fill::randu);
      mat M = C.col_as_mat(0);  // size n x (n + 1)

      writable::list res(5);
      res[0] = as_doubles_matrix(C.slice(0));
      res[1] = as_doubles_matrix(C.slice(1));
      res[2] = as_doubles_matrix(C.slice(2));
      res[3] = as_doubles_matrix(C.slice(3));
      res[4] = as_doubles_matrix(M);

      res.attr("names") = strings({"slice0", "slice1", "slice2", "slice3",
        "col_as_mat"});

      return res;
    }

# 52 Convert column to matrix

`.row_as_mat(row_number)` is a member function of the `Cube` class, it returns a matrix of the specified cube row and the number of columns is preserved. Given a cube of size `R x C x S`, the resultant matrix size is `S x C`.

## 52.1 Examples

    [[cpp4r::register]] list row_as_mat1_(const int& n) {
      cube C(n, n + 1, n + 2, fill::randu);
      mat M = C.row_as_mat(0);  // size (n + 2) x (n + 1)

      writable::list res(5);
      res[0] = as_doubles_matrix(C.slice(0));
      res[1] = as_doubles_matrix(C.slice(1));
      res[2] = as_doubles_matrix(C.slice(2));
      res[3] = as_doubles_matrix(C.slice(3));
      res[4] = as_doubles_matrix(M);

      res.attr("names") = strings({"slice0", "slice1", "slice2", "slice3",
        "row_as_mat"});

      return res;
    }

# 53 Convert sparse matrix to dense matrix

`.as_dense()` is a member function of the `SpMat` class, it avoids the construction of an intermediate sparse matrix representation of the expression.

## 53.1 Examples

    [[cpp4r::register]] doubles as_dense1_(const int& n) {
      sp_mat A;
      A.sprandu(n, n, 0.1);

      // extract column 1 of A directly into dense column vector
      colvec c = A.col(0).as_dense();

      // store the sum of each column of A directly in dense row vector
      rowvec r = sum(A).as_dense();

      return as_doubles(c + r.t());
    }

# 54 Dense matrix and vector transposition

`.t()` is a member function of the `Mat`, `Col` and `Row` classes, it returns a transposed copy of the object. For real matrices, the transpose is a simple transposition of the elements. For complex matrices, the transpose is a Hermitian conjugate transposition of the elements (e.g., the signs of the imaginary components are flipped).

## 54.1 Examples

    [[cpp4r::register]] doubles_matrix<> transpose1_(const int& n) {
      mat A(n, n + 1, fill::randu);
      mat B = A.t();
      return as_doubles_matrix(B);
    }

# 55 Sparse matrix transposition

`.st()` is a member function of the `SpMat` classe, it returns a transposed copy of the object. For real matrices, it is not applicable. For complex matrices, the transpose is a simple transposition of the elements (e.g., the signs of imaginary components are not flipped).

## 55.1 Examples

    [[cpp4r::register]] doubles_matrix<> transpose2_(const int& n) {
      sp_mat A;
      A.sprandu(n, n + 1, 0.1);
      sp_mat B = A.t();
      return as_doubles_matrix(B);
    }

# 56 Matrix inversion

`.i()` is a member function of the `Mat` class, it provides an inverse of the matrix. If the matrix is not square sized, a `std::logic_error` exception is thrown. If the matrix appears to be singular, the output matrix is reset and a `std::runtime_error` exception is thrown.

## 56.1 Examples

    [[cpp4r::register]] doubles inverse1_(const doubles_matrix<>& a,
                                          const doubles b) {
      mat A = as_Mat(a);
      vec B = as_Col(b);

      mat X = inv(A);
      vec Y = X * B;

      return as_doubles(Y);
    }

## 56.2 Caveats

- If the matrix is known to be symmetric positive definite, `inv_sympd()`.
- To solve a system of linear equations, such as `Z = inv(X) * Y`, `solve()` can be faster and/or more accurate.

# 57 Maximum and minimum

`.min()` and `.max()` are member functions of the `Mat`, `Col`, `Row`, and `Cube` classes. These return the minimum and maximum values of the object, respectively. For objects with complex numbers, absolute values are used for comparison.

## 57.1 Examples

    [[cpp4r::register]] doubles maxmin1_(const int& n) {
      mat A = randu<mat>(n, n);

      writable::doubles res(2);
      res[0] = A.max();
      res[1] = A.min();

      res.attr("names") = strings({"max", "min"});

      return res;
    }

# 58 Linear index of maximum and minimum

`.index_min()` and `.index_max()` are member functions of the `Mat`, `Col`, `Row`, and `Cube` classes. They return the linear index of the minimum and maximum values of the object, respectively. For objects with complex numbers, absolute values are used for comparison. The returned index is of type `uword`.

## 58.1 Examples

    [[cpp4r::register]] doubles index_maxmin1_(const int& n) {
      mat A = randu<mat>(n, n);

      writable::doubles res(6);
      res[0] = static_cast<int>(A.index_max());
      res[1] = static_cast<int>(A.index_min());
      res[2] = A(0, 0);
      res[3] = A(1, 0);
      res[4] = A(0, 1);
      res[5] = A(1, 1);

      res.attr("names") = strings({"index_max", "index_min", "element0", "element1",
        "element2", "element3"});

      return res;
    }

# 59 In-range

`.in_range(** i **)` is a member function of `Mat`, `Col`, `Row`, `Cube`, `SpMat` and `field`, it returns `true` if the given location or span is currently valid and `false` if the object is empty, the location is out of bounds, or the span is out of bounds.

Function | Mat | Col | Row | Cube | SpMat | Field
--- | --- | --- | --- | --- | --- | ---
`.in_range(span(start, end))` | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
`.in_range(row, col)` | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
`.in_range(span(start_row, end_row), span(start_col, end_col))` | ✓ | ✓ | ✓ |  | ✓ | ✓
`.in_range(row, col, slice)` |  |  |  | ✓ |  | ✓
`.in_range(span(start_row, end_row), span(start_col, end_col), span(start_slice, end_slice))` |  |  |  | ✓ |  | ✓
`.in_range(first_row, first_col, size(X))` (`X` is a matrix or field) | ✓ | ✓ | ✓ |  | ✓ | ✓
`.in_range(first_row, first_col, size(n_rows, n_cols))` | ✓ | ✓ | ✓ |  | ✓ | ✓
`.in_range(first_row, first_col, first_slice, size(Q))` (`Q` is a cube or field) |  |  |  | ✓ |  | ✓
`.in_range(first_row, first_col, first_slice, size(n_rows, n_cols, n_slices))` |  |  |  | ✓ |  | ✓

Instances of `span(a,b)` can be replaced by:

- `span()` or `span::all` to indicate the entire range.
- `span(a)` to indicate a particular row, column, or slice.

## 59.1 Examples

    [[cpp4r::register]] logicals in_range1_(const int& n) {
      mat A(n, n + 1, fill::randu);

      writable::logicals res(3);
      res[0] = A.in_range(0, 0);
      res[1] = A.in_range(3, 4);
      res[2] = A.in_range(4, 5);

      res.attr("names") = strings({"in_range00", "in_range34", "in_range45"});

      return res;
    }

# 60 Is empty

`.is_empty()` is a member function of the `Mat`, `Col`, `Row`, `Cube`, `SpMat`, and `field` classes. It returns `true` if the object has no elements and `false` if the object has one or more elements.

## 60.1 Examples

    [[cpp4r::register]] logicals is_empty1_(const int& n) {
      mat A(n, n + 1, fill::randu);

      writable::logicals res(2);
      res[0] = A.is_empty();

      A.reset();
      res[1] = A.is_empty();

      res.attr("names") = strings({"before_reset", "after_reset"});

      return res;
    }

# 61 Is vector/column vector/row vector

`.is_vec()`, `.is_colvec()` and `.is_rowvec()` are member functions of `Mat` and `SpMat`.

- `.is_vec()` returns `true` if the matrix can be interpreted as a vector (either column or row vector) and `false` otherwise.
- `.is_colvec()` returns `true` if the matrix can be interpreted as a column vector and `false` otherwise.
- `.is_rowvec()` returns `true` if the matrix can be interpreted as a row vector and `false` otherwise.

## 61.1 Examples

    [[cpp4r::register]] logicals is_vec1_(const int& n) {
      mat A(n, 1, fill::randu);
      mat B(1, n, fill::randu);
      mat C(0, 1, fill::randu);
      mat D(1, 0, fill::randu);

      writable::logicals res(5);
      res[0] = A.is_vec();
      res[1] = A.is_colvec();
      res[2] = B.is_rowvec();
      res[3] = C.is_colvec();
      res[4] = D.is_rowvec();

      res.attr("names") = strings({"Nx1_is_vec", "Nx1_is_colvec", "1xN_is_rowvec",
        "0x1_is_colvec", "1x0_is_rowvec"});

      return res;
    }

## 61.2 Caveat

Do not assume that the vector has elements if these functions return `true`. It is possible to have an empty vector (e.g., 0x1 as in the examples).

# 62 Is sorted

`.is_sorted()`, `.is_sorted(sort_direction)` and `.is_sorted(sort_direction, dim)` are member function of `Mat`, `Row`, and `Col`. For matrices and vectors with complex numbers, order is checked via absolute values.

If the object is a vector, these return a `bool` indicating whether the elements are sorted. If the object is a matrix, these return a `bool` indicating whether the elements are sorted in each column (`dim = 0`, default) or each row (`dim = 1`), and the `dim` argument is optional.

The `sort_direction` argument is optional, `sort_direction` can be one of the following strings:

- `"ascend"`: the elements are ascending, consecutive elements can be equal, and this is the default operation.
- `"descend"`: the elements are descending, and consecutive elements can be equal.
- `"strictascend"`: the elements are strictly ascending, and consecutive elements cannot be equal.
- `"strictdescend"`: the elements are strictly descending, and consecutive elements cannot be equal.

## 62.1 Examples

    [[cpp4r::register]] logicals is_sorted1_(const int& n) {
      vec a(n, fill::randu);
      vec b = sort(a);
      mat A(10, 10, fill::randu);

      writable::logicals res(4);
      res[0] = a.is_sorted();
      res[1] = b.is_sorted();
      res[2] = A.is_sorted("descend", 1);
      res[4] = A.is_sorted("ascend", 1);

      res.attr("names") = strings({"a_sorted", "b_sorted", "A_descend",
        "A_ascend"});

      return res;
    }

# 63 Is upper triangular/lower triangular

`.is_trimatu()` and `.is_trimatl()` are member functions of `Mat` and `SpMat`. `.is_trimatu()` returns `true` if the matrix is upper triangular (e.g., the matrix is square sized and all elements below the main diagonal are zero) and `false` otherwise. `.is_trimatl()` returns `true` if the matrix is lower triangular (e.g., the matrix is square sized and all elements above the main diagonal are zero) and `false` otherwise.

## 63.1 Examples

    [[cpp4r::register]] logicals is_triangular1_(const int& n) {
      mat A(n, n, fill::randu);
      mat B = trimatl(A);

      writable::logicals res(3);
      res[0] = B.is_trimatu();
      res[1] = B.is_trimatl();

      B.reset();
      res[2] = B.is_trimatu();

      res.attr("names") = strings({"is_trimatu", "is_trimatl",
        "is_trimatu_after_reset"});

      return res;
    }

## 63.2 Caveat

If these functions return `true`, do not assume that the matrix contains non-zero elements on or above/below the main diagonal. It is possible to have an empty matrix (e.g., 0x0 as in the examples).

# 64 Is diagonal

`is_diagmat()` is a member function of `Mat` and `SpMat`. It returns `true` if the matrix is diagonal (e.g., all elements outside of the main diagonal are zero). If the matrix is not square sized, a `std::logic_error` exception is thrown.

## 64.1 Examples

    [[cpp4r::register]] logicals is_diagonal1_(const int& n) {
      mat A(n, n, fill::randu);
      mat B = diagmat(A);

      writable::logicals res(3);
      res[0] = A.is_diagmat();
      res[1] = B.is_diagmat();

      A.reset();
      res[2] = A.is_diagmat();

      res.attr("names") = strings({"A_diagmat", "B_diagmat",
        "A_diagmat_after_reset"});

      return res;
    }

## 64.2 Caveat

If this function returns `true`, do not assume that the matrix contains non-zero elements on the main diagonal only. It is possible to have an empty matrix (e.g., 0x0 as in the examples).

# 65 Is square

`.is_square()` is a member function of the `Mat` and `SpMat` classes. It returns `true` if the matrix is square sized (e.g., the number of rows is equal to the number of columns) and `false` otherwise.

## 65.1 Examples

    [[cpp4r::register]] logicals is_square1_(const int& n) {
      mat A(n, n, fill::randu);
      mat B = diagmat(A);

      writable::logicals res(3);
      res[0] = A.is_square();
      res[1] = B.is_square();

      A.reset();
      res[2] = A.is_square();

      res.attr("names") = strings({"A_square", "B_square",
        "A_square_after_reset"});

      return res;
    }

## 65.2 Caveats

If this function returns `true`, do not assume that the matrix is non-empty. It is possible to have an empty matrix (e.g., 0x0 as in the examples).

# 66 Is symmetric

`.is_symmetric()` is a member function of the `Mat` and `SpMat` classes. It returns `true` if the matrix is symmetric (e.g., the matrix is square sized and the transpose is equal to the original matrix) and `false` otherwise.

## 66.1 Examples

    [[cpp4r::register]] logicals is_symmetric1_(const int& n) {
      mat A(n, n, fill::randu);
      mat B = symmatu(A);

      writable::logicals res(3);
      res[0] = A.is_symmetric();
      res[1] = B.is_symmetric();

      A.reset();
      res[2] = A.is_symmetric();

      res.attr("names") = strings({"A_symmetric", "B_symmetric",
        "A_symmetric_after_reset"});

      return res;
    }

## 66.2 Caveats

If this function returns `true`, do not assume that the matrix is non-empty. It is possible to have an empty matrix (e.g., 0x0 as in the examples).

# 67 Is hermitian

`.is_hermitian()` is a member function of the `Mat` and `SpMat` classes. It returns `true` if the matrix is Hermitian or self-adjoint (e.g., the matrix is square sized and the conjugate transpose is equal to the original matrix) and `false` otherwise.

## 67.1 Examples

    [[cpp4r::register]] logicals is_hermitian1_(const int& n) {
      cx_mat A(n, n, fill::randu);
      cx_mat B = A.t() * A;

      writable::logicals res(3);
      res[0] = A.is_hermitian();
      res[1] = B.is_hermitian();

      A.reset();
      res[2] = A.is_hermitian();

      res.attr("names") = strings({"A_hermitian", "B_hermitian",
        "A_hermitian_after_reset"});

      return res;
    }

## 67.2 Caveats

If this function returns `true`, do not assume that the matrix is non-empty. It is possible to have an empty matrix (e.g., 0x0 as in the examples).

# 68 Is symmetric/hermitian positive definite

`.is_sympd()` and `.is_sympd(tol)` are a member function of the `Mat` and `SpMat` classes. It returns `true` if the matrix is symmetric/hermitian positive definite within a tolerance (e.g., the matrix is square sized and all its eigenvalues are positive) and `false` otherwise. The `tol` argument is optional, the default is `tol = 100 * datum::eps * norm(X, "fro")`.

## 68.1 Examples

    [[cpp4r::register]] logicals is_sympd1_(const int& n) {
      mat A(n, n, fill::randu);
      mat B = A * A.t();

      writable::logicals res(3);
      res[0] = A.is_sympd();
      res[1] = B.is_sympd();

      A.reset();
      res[2] = A.is_sympd();

      res.attr("names") = strings({"A_sympd", "B_sympd",
        "A_sympd_after_reset"});

      return res;
    }

# 69 Is zero

`.is_zero()` and `.is_zero(tol)` are a member function of the `Mat`, `Col`, `Row`, `Cube`, and `SpMat` classes. It returns `true` if all elements are zero within a tolerance and `false` otherwise. For complex numbers, each component (real and imaginary) is checked separately. The `tol` argument is optional.

## 69.1 Examples

    [[cpp4r::register]] logicals is_zero1_(const int& n) {
      mat A(n, n, fill::randu);
      cube B(n, n, n, fill::zeros);
      sp_mat C(n, n);

      writable::logicals res(3);
      res[0] = A.is_zero(0.005);
      res[1] = B.is_zero(0.005);
      res[2] = C.is_zero(0.005);

      res.attr("names") = strings({"A_is_zero", "B_is_zero", "C_is_zero"});

      return res;
    }

# 70 Is finite

`.is_finite()` is a member function of the `Mat`, `Col`, `Row`, `Cube`, and `SpMat` classes. It returns `true` if all elements are finite and `false` otherwise.

## 70.1 Examples

    [[cpp4r::register]] logicals is_finite1_(const int& n) {
      mat A(n, n, fill::randu);
      cube B(n, n, n, fill::randu);
      sp_mat C(n, n);

      // Insert infinite values
      B(0, 0, 0) = datum::inf;
      C(0, 0) = -1.0 * datum::inf;

      writable::logicals res(3);
      res[0] = A.is_finite();
      res[1] = B.is_finite();
      res[2] = C.is_finite();

      res.attr("names") = strings({"A_is_finite", "B_is_finite", "C_is_finite"});

      return res;
    }

# 71 Has infinity

`.has_inf()` is a member function of the `Mat`, `Col`, `Row`, `Cube`, and `SpMat` classes. It returns `true` if the object contains at least one infinite value and `false` otherwise.

## 71.1 Examples

    [[cpp4r::register]] logicals has_inf1_(const int& n) {
      mat A(n, n, fill::randu);
      cube B(n, n, n, fill::randu);
      sp_mat C(n, n);

      // Insert infinite values
      B(0, 0, 0) = datum::inf;
      C(0, 0) = -1.0 * datum::inf;

      writable::logicals res(3);
      res[0] = A.has_inf();
      res[1] = B.has_inf();
      res[2] = C.has_inf();

      res.attr("names") = strings({"A_has_inf", "B_has_inf", "C_has_inf"});

      return res;
    }

# 72 Has not-a-number

`.has_nan()` is a member function of the `Mat`, `Col`, `Row`, `Cube`, and `SpMat` classes. It returns `true` if the object contains at least one not-a-number (NaN) value and `false` otherwise.

## 72.1 Examples

    [[cpp4r::register]] logicals has_nan1_(const int& n) {
      mat A(n, n, fill::randu);
      cube B(n, n, n, fill::randu);
      sp_mat C(n, n);

      // Insert NaN values
      B(0, 0, 0) = datum::nan;
      C(0, 0) = -1.0 * datum::nan;

      writable::logicals res(3);
      res[0] = A.has_nan();
      res[1] = B.has_nan();
      res[2] = C.has_nan();

      res.attr("names") = strings({"A_has_nan", "B_has_nan", "C_has_nan"});

      return res;
    }

## 72.2 Caveat

`NaN` is not equal to anything, even itself.