@QuantStack
@QuantStack
Going native: C++ as a first-class citizen of the Jupyter ecosystem
The team
From day 1, Jupyter was developed by scientists for scientists and educators
Programming languages are not only used to build complex systems, but also to explore and gain insight about
Interactive workflows:
(Outside of ROOT) We lack a good story for interative computing in C++
This hurts productivity of C++ software developers.
A C++ kernel for Jupyter based on
A modern C++ implementation of the Jupyter protocol
Cling, the C++ interpreter developed at CERN
Another area where Jupyter shines is Jupyter interactive widgets
Jupyter interactive widgets: thick front-end and thin back-end
Next goal: implement C++ back-end for other widgets libraries
The Lazy Tensor Algebra Expression System
Xtensor is a flexible expression system in which any data structure can be plugged, offering the most expressive API to the users.
A C++ template library for multi-dimensional array manipulation
Followings the idioms of the C++ STL
(iterator pairs, clear value semantics)
Python bindings to enable xtensor APIs on numpy arrays.
Julia bindings to enable xtensor APIs on Julia arrays.
R bindings to enable xtensor APIs on R arrays.
Cookiecutter projects for authoring of Python, Julia, and R extensions.
BLAS bindings to enable BLAS operations on xtensor expressions.
SIMD acceleration kernels.
Ever heard of numpy ?
Python 3 - numpy
C++ 14 - xtensor
np.array([[3, 4], [5, 6]])
arr.reshape([3, 4])
xt::xarray<double>({{3, 4}, {5, 6}})
xt::xtensor<double, 2>({{3, 4}, {5, 6}})
arr.reshape({3, 4});
np.linspace(1.0, 10.0, 100)
np.logspace(1.0, 10.0, 100)
np.arange(3, 7)
np.eye(4)
np.zeros([3, 4])
np.ones([3, 4])
xt::linspace<double>(1.0, 10.0, 100)
xt::logspace<double>(1.0, 10.0, 100)
xt::arange(3, 7)
xt::eye(4)
xt::zeros<double>({3, 4})
xt::ones<double>({3, 4})
Python 3 - numpy
C++ 14 - xtensor
a[:, np.newaxis]
a[:5, 1:]
a[5:1:-1, :]
np.broadcast(a, [4, 5, 7])
np.vectorize(f)
a[a > 5]
a[[0, 1], [0, 0]]
xt::view(a, xt::all(), xt::newaxis())
xt::view(a, xt::range(_, 5), xt::range(1, _))
xt::view(a, xt::range(5, 1, -1), xt::all())
xt::broadcast(a, {4, 5, 7})
xt::vectorize(f)
xt::filter(a, a > 5)
xt::index_view(a, {{0, 0}, {1, 0}})
np.sum(a, axis=[0, 1])
np.sum(a)
np.prod(a, axis=1)
np.prod(a)
np.mean(a, axis=1)
np.mean(a)
xt::sum(a, {0, 1})
xt::sum(a)
xt::prod(a, {1})
xt::prod(a)
xt::mean(a, {1})
xt::mean(a)
Python 3 - numpy
C++ 14 - xtensor
np.where(a > 5, a, b)
np.where(a > 5)
np.any(a)
np.all(a)
np.logical_and(a, b)
np.logical_or(a, b)
xt::where(a > 5, a, b)
xt::where(a > 5)
xt::any(a)
xt::all(a)
a && b
a || b
np.absolute(a)
np.exp(a)
np.sqrt(a)
np.cos(a)
np.cosh(a)
scipy.special.erf(a)
np.isnan(a)
xt::abs(a)
xt::exp(a)
xt::sqrt(a)
xt::cos(a)
xt::cosh(a)
xt::erf(a)
xt::isnan(a)
Python 3 - numpy
C++ 14 - xtensor
np.random.seed(0)
np.random.randn(10, 10)
np.random.randint(10, 10)
np.random.rand(3, 4)
xt::random::seed(0)
xt::random::randn<double>({10, 10})
xt::random::randint<int>({10, 10}})
xt::random::rand<double>({3, 4}})
np.stack([a, b, c], axis=1)
np.concatenate([a, b, c], axis=1)
xt::stack(xtuple(a, b, c), 1)
xt::concatenate(xtuple(a, b, c), 1)
Broadcasting
xarray<int> a = {{1, 2, 3, 4},
{5, 6, 7, 8},
{9, 10, 11, 12}};
xarray<int> b = { 1, 3, 5, 7};
xarray<int> res = a + b;
Broadcasting
xarray<double> a = {1., 2., 3., 4.};
auto res = xt::broadcast(a, {3, 4});
xarray<double> a, b, c;
// ... initialization of a, b, and c ...
auto res = broadcast(a + b * c, {3, 4});
Iteration
for x in np.nditer(a)
for(auto it=a.begin(); it!=a.end(); ++it)
a with a prescribed broadcasting shape
a.begin({3, 4})
a.end({3, 4})
a in a column-major fashion
a.template begin<layout_type::column_major>()
a.template end<layout_type::column_major>()
a in a column-major fashion with a prescribed broadcasting shape
a.template begin<layout_type::column_major>({3, 4})
a.template end<layout_type::column_major>({3, 4})
http://quantstack.net/xtensor
A Simple Python extension (1/2)
C++: Using an algorithm from the STL on a numpy array
#include <numeric> // Standard library import for std::accumulate
#include "pybind11/pybind11.h" // Pybind11 import to define Python bindings
#include "xtensor/xmath.hpp" // xtensor import for the C++ universal functions
#define FORCE_IMPORT_ARRAY // numpy C api loading
#include "xtensor-python/pyarray.hpp" // Numpy bindings
double sum_of_sines(xt::pyarray<double>& m)
{
auto sines = xt::sin(m);
// sines does not actually hold any value, which are only computed upon access
return std::accumulate(sines.begin(), sines.end(), 0.0);
}
PYBIND11_PLUGIN(xtensor_python_test)
{
xt::import_numpy();
pybind11::module m("xtensor_python_test", "Test module for xtensor python bindings");
m.def("sum_of_sines", sum_of_sines,
"Computes the sum of the sines of the values of the input array");
return m.ptr();
}
A simple Python extension (1/2)
Python: Using an algorithm from the STL on a numpy array
import numpy as np
import xtensor_python_test as xt
a = np.arange(15).reshape(3, 5)
xt.sum_of_sines(v)
A simple Python extension (2/2)
C++: Create a universal function from a C++ scalar function
#include "pybind11/pybind11.h"
#define FORCE_IMPORT_ARRAY
#include "xtensor-python/pyvectorize.hpp"
#include <numeric>
#include <cmath>
namespace py = pybind11;
double scalar_func(double i, double j)
{
return std::sin(i) - std::cos(j);
}
PYBIND11_PLUGIN(xtensor_python_test)
{
xt::import_numpy();
py::module m("xtensor_python_test", "Test module for xtensor python bindings");
m.def("vectorized_func", xt::pyvectorize(scalar_func), "");
return m.ptr();
}
A simple Python extension (2/2)
Python: Create a numpy-style universal function from a C++ scalar function
import numpy as np
import xtensor_python_test as xt
x = np.arange(15).reshape(3, 5)
y = [1, 2, 3, 4, 5]
xt.vectorized_func(x, y)
A Simple Julia extension (1/2)
C++: Using an algorithm from the STL on a Julia array
#include <numeric> // Standard library import for std::accumulate
#include "jlcxx/jlcxx.hpp¨ // CxxWrap import to define Julia bindings
#include "xtensor-julia/jltensor.hpp" // Import the jltensor container definition
#include "xtensor/xmath.hpp" // xtensor import for the C++ universal functions
double sum_of_sines(xt::jltensor<double, 2> m)
{
auto sines = xt::sin(m); // sines does not actually hold values.
return std::accumulate(sines.cbegin(), sines.cend(), 0.0);
}
JULIA_CPP_MODULE_BEGIN(registry)
jlcxx::Module mod = registry.create_module("xtensor_julia_test");
mod.method("sum_of_sines", sum_of_sines);
JULIA_CPP_MODULE_END
A simple Julia extension (1/2)
Julia: Using an algorithm from the STL on a Julia array
using xtensor_julia_test
arr = [[1.0 2.0]
[3.0 4.0]]
sum_of_sines(arr)
A simple Julia extension (2/2)
C++: Create a numpy-style universal function from a C++ scalar function
#include "jlcxx/jlcxx.hpp"
#include "xtensor-julia/jlvectorize.hpp"
double scalar_func(double i, double j)
{
return std::sin(i) - std::cos(j);
}
JULIA_CPP_MODULE_BEGIN(registry)
jlcxx::Module mod = registry.create_module("xtensor_julia_test");
mod.method("vectorized_func", xt::jlvectorize(scalar_func));
JULIA_CPP_MODULE_END
A simple Julia extension (2/2)
Julia: Create a numpy-style universal function from a C++ scalar function
using xtensor_julia_test
x = [[ 0.0 1.0 2.0 3.0 4.0]
[ 5.0 6.0 7.0 8.0 9.0]
[10.0 11.0 12.0 13.0 14.0]]
y = [1.0, 2.0, 3.0, 4.0, 5.0]
xt.vectorized_func(x, y)
A Simple R extension
C++: Using an algorithm from the STL on a R array
// [[Rcpp::export]]
double sum_of_sines(xt::rarray<int> m)
{
auto sines = xt::sin(m); // sines does not actually hold values.
return std::accumulate(sines.cbegin(), sines.cend(), 0.0);
}
A simple R extension
R: Using an algorithm from the STL on a R array
library('xtensor_r_test')
arr <- array(c(c(1, 2), c(3, 4)), dim = c(2, 2)
xtensor_r_test::sum_of_sines(arr)
Generate your packaged xtensor extension
Bindings with BLAS libraries
BLAS-based implementation of numpy.linalg
numpy.linalg is main goalWe benchmark.
A lot!
GitHub
Documentation
