Examples

On a high level CARMA provides two ways to work with Numpy arrays and Armadillo, see the Function specifications section for details about the available functions and the examples directory for runnable examples.

Manual conversion

The easiest way to use CARMA is manual conversion, it gives you the most control over when to copy or not. You pass a Numpy array as an input and/or as the return type and call the respective conversion function.

Warning

Carma will avoid copying by default so make sure not to return the memory of the input array without copying. If you don’t copy out, the memory is aliased by both the input and output arrays which will cause a segfault.

Borrow

#include <carma>
#include <armadillo>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

py::array_t<double> manual_example(py::array_t<double> & arr) {
    // convert to armadillo matrix without copying.
    // Note the size of the matrix cannot be changed when borrowing
    arma::Mat<double> mat = carma::arr_to_mat<double>(arr);

    // normally you do something useful here ...
    arma::Mat<double> result = arma::Mat<double>(arr.shape(0), arr.shape(1), arma::fill::randu);

    // convert to Numpy array and copy out
    return carma::mat_to_arr(result, true);
}

Transfer ownership

If you want to transfer ownership to the C++ side you can use:

#include <carma>
#include <armadillo>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

py::array_t<double> steal_array(py::array_t<double> & arr) {
    // convert to armadillo matrix
    arma::Mat<double> mat = carma::arr_to_mat<double>(std::move(arr));
    // armadillo now owns and manages the lifetime of the memory
    // if you want to return you don't need to copy out
    return mat_to_arr(mat);
}

py::array_t<double> copy_array(py::array_t<double> & arr) {
    // convert to armadillo matrix
    arma::Mat<double> mat = carma::arr_to_mat<double>(arr, true);
    // armadillo now owns and manages the lifetime of the memory
    // if you want to return you don't need to copy out
    return mat_to_arr(mat);
}

py::array_t<double> copy_const_array(const py::array_t<double> & arr) {
    // convert to armadillo matrix
    arma::Mat<double> mat = carma::arr_to_mat<double>(arr);
    // armadillo now owns and manages the lifetime of the memory
    // if you want to return you don't need to copy out
    return mat_to_arr(mat);
}

Automatic conversion

For automatic conversion you specify the desired Armadillo type for either or both the return type and the function parameter. When calling the function from Python, Pybind11 will call CARMA’s type caster when a Numpy array is passed or returned, see Return policies for details.

Warning

Make sure to include carma in every compilation unit that makes use of the type caster, not including it results in undefined behaviour.

#include <carma>
#include <armadillo>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

arma::Mat<double> automatic_example(arma::Mat<double> & mat) {
    // normally you do something useful here with mat ...
    arma::Mat<double> rand = arma::Mat<double>(mat.n_rows, mat.n_cols, arma::fill::randu);

    arma::Mat<double> result = mat + rand;
    // type caster will take care of casting `result` to a Numpy array.
    return result;
}

Warning

The automatic conversion will not copy the Numpy array’s memory when converting to Armadillo objects. When converting back to Numpy arrays the memory will not be copied out by default. You shoud specify return_value_policy::copy if you want to return the input array.

ArrayStore

There are use-cases where you would want to keep the data in C++ and only return when requested. For example, you write an Ordinary Least Squares (OLS) class and you want to store the residuals, covariance matrix, … in C++ for when additional tests need to be run on the values without converting back and forth.

ArrayStore is a convenience class that provides conversion methods back and forth. It is intended to be used as an attribute as below:

Warning

The ArrayStore owns the data, the returned numpy arrays are views that are tied to the lifetime of ArrayStore.

#include <armadillo>
#include <carma>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

class ExampleClass {
    private:
        carma::ArrayStore<arma::Mat<double>> _x;
        carma::ArrayStore<arma::Mat<double>> _y;

    public:
        ExampleClass(py::array_t<double> & x, py::array_t<double> & y) :
        // copy the array and store it as an Armadillo matrix
        _x{carma::ArrayStore<arma::Mat<double>>(x, true)},
        // steal the array and store it as an Armadillo matrix
        _y{carma::ArrayStore<arma::Mat<double>>(y, false)},

        py::array_t<double> member_func() {
            // normallly you would something useful here
            _x.mat += _y.mat;
            // return mutable view off arma matrix
            return _x.get_view(true);
        }
};

void bind_exampleclass(py::module &m) {
    py::class_<ExampleClass>(m, "ExampleClass")
        .def(py::init<py::array_t<double> &, py::array_t<double> &>(), R"pbdoc(
            Initialise ExampleClass.

            Parameters
            ----------
            arr1: np.ndarray
                array to be stored in armadillo matrix
            arr2: np.ndarray
                array to be stored in armadillo matrix
        )pbdoc")
        .def("member_func", &ExampleClass::member_func, R"pbdoc(
            Compute ....
        )pbdoc");
}

Ordinary Least Squares

Combining the above approaches to compute the Ordinary Least Squares:

#include <carma>
#include <armadillo>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/pytypes.h>

py::tuple ols(arma::mat& X, arma::colvec& y) {
    // We borrow the data underlying the numpy arrays
    int n = X.n_rows, k = X.n_cols;

    arma::colvec coeffs = arma::solve(X, y);
    arma::colvec resid = y - X * coeffs;

    double sig2 = arma::as_scalar(arma::trans(resid) * resid / (n-k));
    arma::colvec std_errs = arma::sqrt(sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );

    // We take ownership of the memory from the armadillo objects and
    // return to python as a tuple containing two Numpy arrays.
    return py::make_tuple(
        carma::col_to_arr(coeffs),
        carma::col_to_arr(std_errs)
    );
}

// adapted from https://gallery.rcpp.org/articles/fast-linear-model-with-armadillo/

Which can be called using:

y = np.linspace(1, 100, num=100) + np.random.normal(0, 0.5, 100)
X = np.hstack((np.ones(100)[:, None], np.arange(1, 101)[:, None]))
coeff, std_err = carma.ols(X, y)

The repository contains tests, examples and CMake build instructions that can be used as an reference.