Expressions¶
a * 2 + b does not compute anything. It builds an expression: a small object that remembers the operands and the operations, and reads values when somebody asks for them. Making a tensor of it, assigning it, reading one element or summing it are the moments something is computed, and then it is computed once, over memory, with no temporary tensors in between.
Every code block on this page is part of examples/13_expressions.cpp, which is built and run with the tests. The output shown is its output.
using strata::Dyn;
using strata::Matrix;
using strata::MatrixX;
using strata::Shape;
using strata::Tensor;
Lazy by default¶
const Matrix<float, 2, 3> a =
strata::generate(Shape<2, 3>{}, [](std::size_t i, std::size_t j) { return static_cast<float>(i + (2 * j)); });
const MatrixX<float> b{Shape<Dyn, Dyn>{2, 3}, 0.5F};
// An expression records what to compute, and computes nothing.
const auto e = (a * 2.0F) + b;
// A tensor is made of it: computed now, in one pass over memory. The static shape wins.
Matrix c = e;
static_assert(std::same_as<decltype(c), Tensor<float, Shape<2, 3>>>);
std::println("c:\n{::.1f}", c);
// Assigned: the extents agree, so it is written in place.
c = c - b;
std::println("c - b:\n{::.1f}", c);
Three things happened. e was built: no loop ran. Matrix c = e ran one loop and wrote the result into new memory, without filling it with zeros first. And c = c - b found that c already has the extents of the result and wrote into it.
The type of the result follows the expression: here float values and the shape both operands agree on, static wherever either fixes it. An expression borrows a tensor that is named, so it never copies one; a temporary tensor, whose elements would be gone before the expression is computed, is moved into the expression, which then owns it.
Shapes are checked. Static extents that differ do not compile; dynamic ones are compared when the expression is made.
template <typename L, typename R>
concept addable = requires(const L& lhs, const R& rhs) { lhs + rhs; };
static_assert(addable<Tensor<float, Shape<2, 3>>, Tensor<float, Shape<Dyn, Dyn>>>);
static_assert(!addable<Tensor<float, Shape<2, 3>>, Tensor<float, Shape<3, 2>>>); // static extents differ
Masks and where¶
A comparison gives a mask: an expression of bool. Masks combine element by element with &, |, ^ and !. The operators && and || are left as they are, because an overloaded one would no longer short-circuit. == stays a single bool, as std::equality_comparable expects; strata::equal is the element-wise one.
// A mask: a bool for every element.
const auto big = a > 2.0F;
// 2 where the mask is true, else the element of a.
const auto capped = strata::where(big, 2.0F, a);
std::println("big: {:l}", big);
std::println("capped: {:l}", capped);
std::println("{} of {} are big, any: {}", strata::count(big), big.size(), strata::any(big));
big: [[false, false, true], [false, true, true]]
capped: [[0, 2, 2], [1, 2, 2]]
3 of 6 are big, any: true
where(mask, a, b) chooses element by element; either branch may be a tensor or a plain number. count, any and all answer questions about a mask over all its elements at once.
Pipes¶
Every function of strata takes its operand from the left of a |, so a computation reads in the order it happens:
const Tensor<double, Shape<Dyn>> x = strata::linspace(-2.0, 2.0, 5);
// sqrt(abs(x)), read in the order things happen.
std::println("{::.3f}", x | strata::abs | strata::sqrt);
// `pow(2.0)` waits for its first argument.
std::println("{::.3f}", x | strata::pow(2.0));
// A pipeline is a value of its own.
constexpr auto norm = strata::square | strata::sum | strata::sqrt;
std::println("norm(x) = {:.3f}", x | norm);
A call that stops short of its first argument keeps the rest and waits: pow(2.0) is pow(x, 2.0) once x arrives. That makes a pipeline a value of its own, norm above, to apply to any tensor or pass to a function.
| binds looser than the comparisons, so a comparison after a pipe needs parentheses, (x | abs) > 1. And between two masks | is still the bitwise or: only a function on its right makes a pipe.
Reductions¶
A reduction turns many values into one. Of a whole expression it gives a value at once; along axes it gives an expression of the axes that are left.
const Tensor<double, Shape<2, 3>> m =
strata::generate(Shape<2, 3>{}, [](std::size_t i, std::size_t j) { return static_cast<double>(i + (2 * j)); });
// Of the whole tensor: a value, at once.
std::println("sum {}, mean {}, max {}", strata::sum(m), strata::mean(m), strata::max(m));
// Along an axis: an expression of the axes that are left.
const Tensor<double, Shape<2>> per_row = strata::sum.along<1>(m);
const Tensor<double, Shape<3>> per_column = strata::sum.along<0>(m);
std::println("per row {:l}, per column {:l}", per_row, per_column);
The functions are sum, prod, min, max, mean, all, any and count. A sum is made in the type T + T gives, which is why bytes add up without wrapping; sum.as<double>(x) chooses the accumulator, and a float sum is pairwise, so a million small terms do not drift. They are expressions like the rest: slice one, use it as an operand, put it in a pipeline, square | mean.along<1, 2> | sqrt.
Generators¶
Some tensors have values but no memory: a constant, the identity, a run of numbers, a function of the index. They mix into expressions like views and become tensors when assigned; nothing is stored until then.
std::println("iota: {:l}", strata::iota<int>(5));
std::println("linspace: {:l:.2f}", strata::linspace(0.0, 1.0, 5));
std::println("eye:\n{}", strata::eye<int>(Shape<3, 3>{}));
// No memory until the tensor is made.
const Tensor<int, Shape<3, 3>> scaled = strata::eye<int>(Shape<3, 3>{}) * 7;
std::println("7 I: {:l}", scaled);
iota: [0, 1, 2, 3, 4]
linspace: [0.00, 0.25, 0.50, 0.75, 1.00]
eye:
[[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]
7 I: [[7, 0, 0], [0, 7, 0], [0, 0, 7]]
There are also zeros, ones, full, their _like forms that copy the shape and type of a tensor at hand, and generate(shape, f) for any function of the indices. Slicing a generator gives a generator.
einsum¶
einsum writes the products of linear algebra in one notation. Each operand names its axes with letters; the letters after -> are the result's, and every other letter is summed over.
// [[1, 2], [3, 4]].
const Matrix<double, 2, 2> p = strata::generate(
Shape<2, 2>{}, [](std::size_t i, std::size_t j) { return static_cast<double>((2 * i) + j + 1); });
// The matrix product.
std::println("p p:\n{}", Tensor{strata::einsum<"ij,jk->ik">(p, p)});
// No letters left after the arrow: rank 0, read with `[]`.
std::println("trace: {}", strata::einsum<"ii->">(p)[]);
"ij,jk->ik" is the matrix product, "ii->" the trace (a letter twice in one operand takes the diagonal), "i,j->ij" the outer product. The result is an expression like any other: reading one element computes one sum. The einsum guide goes through more of them.
Assignment and aliasing¶
When an expression is assigned to a tensor that it reads, the order matters. t = t * 2.0 reads each element just where it writes it, which is safe. t = t | permute<1, 0> does not: by the time the second element is read, the first write has already overwritten it.
strata works out which it is, mostly at compile time and otherwise by a few pointer comparisons, and for the second kind computes into a buffer first. Nothing in your code says so:
Matrix<double, 2, 2> t = p;
// Reads each element where it writes it: computed in place.
t = t * 2.0;
// Would read elements already written: computed into a buffer first.
t = t | strata::permute<1, 0>;
std::println("transposed, doubled: {:l}", t);
// A promise that nothing is read after it is written: no check, no buffer.
t = t * 2.0 | strata::noalias;
std::println("doubled again: {:l}", t);
If you know better, noalias says so: the check and the buffer are skipped. It is a promise. Where precondition checks are active it is verified; where they are not, a false one is undefined behaviour.
The same goes for the parts of an expression that are read many times and dear to compute again, for instance the operand of an einsum that is itself an expression. The assignment computes them into a buffer first, and the cost model decides, at compile time. scratch chooses where those buffers take their memory from; 11_scratch shows it.
See also¶
- Element-wise expressions, Reductions, Generators, einsum and Evaluation and assignment in the API reference.
- Tensors and views: the things expressions are made of.