-
Notifications
You must be signed in to change notification settings - Fork 5
Home
ra-ra (version 0.1, updated 2016 September 29)
(c) Daniel Llorens 2016
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
ra-ra is a general purpose multidimensional array and expression template library for C++14/C++17. Please remember that this manual is a work in progress. There are many errors and whole sections yet to be written.
| • Overview: | Array programming and C++. | |
| • Usage: | Everything you can do with ra-ra. | |
| • Reference: | Systematic list of types and functions. | |
| • Hazards: | User beware. | |
| • The future: | Could be even better. |
A multidimensional array is a container whose elements can be looked up using a multi-index (i₀, i₁, ...). Each of the indices i₀, i₁, ... has a fixed range [0, n₀), [0, n₁), ... so the array is ‘rectangular’. The number of indices in the multi-index is the rank of the array, and the list (n₀, n₁, ... nᵣ₋₁) is the shape of the array. We speak of a rank-r array or of an r-array.
Often we deal with multidimensional expressions where the elements aren’t stored anywhere, but are computed on demand when the expression is looked up. In this general sense, an ‘array’ is just a function of integers with a rectangular domain.
Arrays (also called matrices, vectors, or tensors) are everywhere in math and many other fields, and it is enormously useful to be able to manipulate arrays as individual entities rather than as aggregates. Not only is
A = B+C;
much more compact and easier to read than
for (int i=0; i!=m; ++i)
for (int j=0; j!=n; ++j)
for (int k=0; k!=p; ++k)
A(i, j, k) = B(i, j, k)+C(i, j, k);
but it’s also safer and less redundant. For example, the order of the loops may be something you don’t really care about.
However, if array operations are implemented naively, a piece of code such as A=B+C may result in the creation of a temporary to hold B+C which is then assigned to A. Needless to say this is very wasteful if the arrays involved are large.
Fortunately the problem is almost as old as aggregate data types, and other programming languages have addressed it with optimizations such as ‘loop fusion’ REF, ‘drag along’ REF, or ‘deforestation’ REF. In the C++ context the technique of ‘expression templates’ was pioneered in the late 90s by libraries such as Blitz++ REF. It works by making B+C into an ‘expression object’ which holds references to its arguments and executes the sum only when its elements are looked up. The compiler removes the temporary expression objects during optimization, so that A=B+C results (in principle) in the same generated code as the complicated loop nest above.
| • Rank polymorphism: | What makes arrays special. | |
| • Drag along and beating: | The basic array optimizations. | |
| • Why C++: | High level, low level. | |
| • Guidelines: | How ra-ra tries to do things. | |
| • Other libraries: | Inspiration and desperation. |
Next: Drag along and beating, Up: Overview
Rank polymorphism is the ability to treat an array of rank r as an array of lower rank where the elements are themselves arrays.
For example, think of a matrix A, a 2-array with sizes (n₀, n₁) where the elements A(i₀, i₁) are numbers. If we consider the subarrays (rows) A(0, ...), A(1, ...), ..., A(n₀-1, ...) as individual elements, then we have a new view of A as a 1-array of size n₀ with those rows as elements. We say that the rows A(i₀)≡A(i₀, ...) are the 1-cells of A, and the numbers A(i₀, i₁) are 0-cells of A. For an array of arbitrary rank r the (r-1)-cells of A are called its items.
An obvious way to store an array in linearly addressed memory is to place its items one after another. So we would store a 3-array as
A: [A(0), A(1), ...]
and the items of A(i₀), etc. are in turn stored in the same way, so
A: [A(0): [A(0, 0), A(0, 1) ...], ...]
and the same for the items of A(i₀, i₁), etc.
A: A(0, 0): [A(0, 0, 0), A(0, 0, 1) ...], A(0, 1): [A(0, 1, 0), A(0, 1, 1) ..., ...]
This way to lay out an array in memory is called row-major order or C-order, since it’s the default order for arrays in C (see Other libraries). A row-major array A with sizes (n₀, n₁, ... nᵣ₋₁) can be looked up like this:
A(i₀, i₁, ...) = (storage-of-A) [(((i₀n₁ + i₁)n₂ + i₂)n₃ + ...)+iᵣ₋₁] = (storage-of-A) [o + s₀i₀ + s₁i₁ + ...]
where the numbers (s₀, s₁, ...) are called the strides or the dope vector REF. Note that the ‘linear’ or ‘raveled’ address [o + s₀i₀ + s₁i₁ + ...] is an affine function of (i₀, i₁, ...). Therefore if we represent an array as a tuple
A ≡ ((storage-of-A), o, (s₀, s₁, ...))
then any affine transformation of the indices can be achieved simply by modifying the numbers (o, (s₀, s₁, ...)), with no need to touch the storage. This includes very common operations such as: transposing axes, reversing the order along an axis, most cases of slicing, and some important cases of reshaping or tiling the array.
A trivial example is the i₀-th item of A:
A(i₀) ≡ ((storage-of-A), o+s₀i₀, (s₁, ...))
Note that we can iterate over these items by simply bumping the pointer o+s₀i₀. This means that iterating over (k>0)-cells doesn’t cost any more than iterating over 0-cells.
Next: Why C++, Previous: Rank polymorphism, Up: Overview
Next: Guidelines, Previous: Drag along and beating, Up: Overview
Of course the main reason is that (this being a personal project) I’m more familiar with C++ than with other languages to which the following applies.
However, C++ is unique in that it supports the low level control that is necessary for interoperation with external libraries and languages, but still has the abstraction power to create the features we want even though the language has no native support for most of them.
The classic array languages, APL and J, have their support for arrays baked in. The same is true for other languages with array facilities such as Fortran or Octave/Matlab. Array libraries for general purpose languages usually depend heavily on C extensions. In Numpy’s case this is both for reasons of flexibility (e.g. to obtain predictable memory layout and machine types) and of performance.
On the other extreme, an array library for C would be hampered by the limited means of abstraction in the language (no polymorphism, no metaprogramming, etc.) so that the natural choice of C programmers is to resort to code generators, which eventually turn into new languages.
In C++, a library is enough.
Next: Other libraries, Previous: Why C++, Up: Overview
ra-ra attempts to be general, consistent, and transparent.
Generality is achieved by removing arbitrary restrictions and by adopting the rank extension mechanism of J. ra-ra supports array operations with an arbitrary number of arguments. Any of the arguments in an array expression can be read from or written to. Arrays or array expressions can be of any rank. Slicing operations work for subscripts of any rank, as in APL. You can freely use your own types as array elements.
Consistency is achieved by having a clear set of concepts and having the realizations of those concepts adhere to the concept as closely as possible. ra-ra offers a few different types of views and containers, but it should be possible to use them interchangeably whenever the properties that justify their existence are not involved. When this isn’t possible, it’s a bug. For example, you can currently create a higher rank iterator on a View but not a SmallView; this is a bug.
Sometimes consistency requires a choice. For example, given array views A and B, A=B copies the contents of view B into view A. To change view A instead (to treat A as a pointer) would be the default meaning of A=B for C++ types, and result in better consistency with the rest of the language, but I have decided that having consistency between views and containers (which ‘are’ their contents in a sense that views aren’t) is more important.
Transparency is achieved by avoiding opaque types. An array view consists of a pointer and a list of strides and I see no point in hiding that. Manipulating the strides directly is often useful. A container consists of storage and a view and that isn’t hidden either. Some of the types have an obscure implementation but I consider that a defect. Ideally you should be able to rewrite expressions on the fly, or plug in your own traversal methods or storage handling.
That isn’t to mean that you need to be in command of a lot of internal detail to be able to use the library. I hope to have provided a high level interface to most operations and a reasonably sweet syntax. However, transparency is critical to achieve interoperation with external libraries and languages. When you need to, you’ll be able to guarantee that your array is stored by compact columns or that the real parts are interleaved with the imaginary parts.
Previous: Guidelines, Up: Overview
Here I try to list the C++ array libraries that I know of, or libraries that I think deserve a mention for the way they deal with arrays. It is not an extensive review, since I have only used a few of these libraries myself. Please follow the links if you want to be properly informed.
Since the C++ standard library doesn’t offer a standard multidimensional array type, some libraries for specific tasks (linear algebra operations, finite elements, optimization) offer an accessory array library, which may be more or less general. Other libraries have generic array interfaces without needing to provide an array type. FFTW is a good example, maybe because it isn’t C++!
C++ offers multidimensional arrays as a legacy feature from C, e.g. int a[3][4]. These decay to pointers when you do nearly anything with them, don’t know their own sizes or rank, and are generally too limited.
The C++ standard library also offers a number of containers that can be used as 1-arrays, of which the most important are <array>, <vector> and <valarray>. Neither supports higher ranks out of the box, but <valarray> offers array operations for 1-arrays. ra-ra makes use of <array> and <vector> for storage and bootstrapping so we’ll mention these containers from time to time.
Blitz++ (REF) pioneered the use of expression templates in C++. It supported higher rank arrays, as high as it was practical in C++ at the time, but not dynamic rank. It also supported small arrays with compile time sizes, and convenience features such as Fortran-order constructors and arbitrary lower bounds for the array indices (both of which ra-ra doesn’t support on purpose). Storage for large arrays was reference-counted, while in ra-ra that is optional but not the default. It placed a strong emphasis on performance, with array traversal methods such as blocking, space filling curves, etc. To date it remains, I believe, one of the most general array libraries for C++. However, the implementation had to fight the limitations of C++98 and it offered no general rank extension mechanism.
TODO More libraries.
TODO Maybe review other languages, at least the important ones (Fortran/APL/J/Matlab/Numpy).
This is an extended exposition of the features of ra-ra and is probably best read in order. For details on specific functions or types, please see Reference.
| • Using ra-ra: | ra-ra is a header-only library. | |
| • Containers and views: | Data objects. | |
| • Array operations: | Building and executing expressions. | |
| • Rank extension: | How array operands are matched. | |
| • Subarray iteration: | At any rank. | |
| • Slicing: | Subscripting is a special operation. | |
| • Special objects: | Not arrays, yet arrays. | |
| • The rank conjunction: | J comes to C++. | |
| • Compatibility: | With the STL and others. | |
| • Extension: | Using your own types and more. | |
| • Internals: | For all the world to see. |
Next: Containers and views, Up: Usage
ra-ra is a header only library with no dependencies, so you just need to place the ra/ folder somewhere in your include path and add #include <ra/operators.H> and <ra/io.H> at the top of your sources.
A C++14 compiler with partial C++17 support is required. At the time of writing this means gcc 6.2 with -std=c++1z. Most tests pass under clang 3.9 with a couple of extra flags (-Wno-missing-braces, -DRA_OPTIMIZE_SMALLVECTOR=0).
Here is a minimal program:
#include <ra/operators.H> #include <ra/io.H> #include <iostream>int main() { ra::Owned<char, 2> A({2, 5}, "helloworld"); std::cout << format_array(transpose<1, 0>(A), false, "|") << std::endl; }
-| h|w e|o l|r l|l d|d
Next: Array operations, Previous: Using ra-ra, Up: Usage
ra-ra offers two kinds of data objects. The first is the container, whose defining characteristic is that it owns its data. Creating a container uses memory and destroying it causes that memory to be freed.
There are three kinds of containers: fixed size, fixed rank/variable size, and variable rank. Here fixed means ‘compile time constant’ while variable is normally a run time constant. (Some variable size arrays can be resized but variable rank arrays cannot normally have their rank changed. Instead, you create a new container or view with the rank you want.)
For example:
{
ra::Small<double, 2, 3> a(0.); // a fixed size 2x3 array
ra::Owned<double, 2> b({2, 3}, 0.); // a variable size 2x3 array
ra::Owned<double> c({2, 3}, 0.); // a variable rank 2x3 array
// a, b, c destroyed at end of scope
}
The main reason to have all these different types is performance; the compiler can do a much better job when it knows the size or the rank of the array. Also, the sizes of a fixed size array do not need to be stored in memory so if you have thousands of small arrays it certainly pays off to use the fixed size types for them. Fixed size or fixed rank arrays are also safer to use; sometimes ra-ra will be able to detect errors in the sizes or ranks of array operands at compile time, if the appropriate types are used.
Container constructors come in two main forms. The first takes a single argument which is copied into the new container. This argument provides shape information if the container type requires it.
ra::Small<double, 2, 3> a = 0.; // 0. is copied into a ra::Small<double, 2, 3> b = a; // the contents of a are copied into b ra::Owned<double> c = a; // c takes the size of a and a is copied into c ra::Owned<double> d = 0.; // d is a 0-array with one element d()==0.
The second form takes two arguments, one giving the shape, the second the contents.
ra::Owned<double, 2> a({2, 3}, 1.); // a has size 2x3 and be filled with 1.
ra::Owned<double> b({2, 3}, a); // b has size 2x3 and a is copied into b
The last example may result in an error if the shape of a and {2, 3} don’t match. Here the shape of 1. (which is ()) matches (2, 3) by a mechanism of rank extension (see Rank extension).
Finally, there are special constructors where the content argument is either a pointer or an std::initializer_list. This argument isn’t used for shape 1, but only as the (row-major) ravel of the content. The pointer constructor is unsafe —use at your own risk! Nested std::initializer_list isn’t supported yet, but may be supported in the future.
ra::Owned<double, 2> a({2, 3}, {1, 2, 3, 4, 5, 6}); // {{1, 2, 3}, {4, 5, 6}}
double bx[6] = {1, 2, 3, 4, 5, 6}
ra::Owned<double, 2> b({3, 2}, bx); // {{1, 2}, {3, 4}, {5, 6}}
double cx[4] = {1, 2, 3, 4}
ra::Owned<double, 2> c({3, 2}, cx); // *** WHO NOSE ***
using sizes = mp::int_list<2, 3>;
using strides = mp::int_list<1, 2>;
ra::SmallArray<real, sizes, strides> a { 1, 2, 3, 4, 5, 6 }; // {{1, 2, 3}, {4, 5, 6}} stored column-major
Sometimes the pointer constructor gets in the way (TODO):
ra::Owned<char const *, 1> A({3}, "hello"); // ERROR: try to convert char to char const *
ra::Owned<char const *, 1> A({3}, ra::scalar("hello")); // ok, "hello" is a single item
std::cout << format_array(A, false, "|") << std::endl;
-| hello|hello|hello
A view is similar to a container in that it points to actual data in memory. However, the view doesn’t own that data and destroying the view won’t affect it. For example:
ra::Owned<double> c({2, 3}, 0.); // a variable rank 2x3 array
{
auto c1 = c(1); // the first row of array c
// c1 is destroyed here
}
// can still use c here
The data accessed through a view is the data of the ‘root’ container, so modifying the first will be reflected in the latter.
ra::Owned<double> c({2, 3}, 0.);
auto c1 = c(1);
c1(2) = 9.; // c(1, 2) = 9.
Just as for containers, there are separate types of views depending on whether the size is known at compile time, the rank is known at compile time but the size is not, or neither the size nor the rank are known at compile time. ra-ra has functions to create the most common kinds of views:
ra::Owned<double> c({2, 3}, {1, 2, 3, 4, 5, 6});
auto ct = transpose<1, 0>(c); // {{1, 4}, {2, 5}, {3, 6}}
auto cr = reverse(c, 0); // {{4, 5, 6}, {1, 2, 3}}
However, views can point to anywhere in memory and that memory doesn’t have to belong to a ra-ra container. For example:
int raw[6] = {1, 2, 3, 4, 5, 6};
ra::View<int> v1({{2, 3}, {3, 1}}, raw); // view with sizes {2, 3} strides {3, 1}
ra::View<int> v2({2, 3}, raw); // same, default C (row-major) strides
Containers can be treated as views and the container types convert implicitly to view types of the same ‘fixedness’. Therefore if you declare a function
void f(ra::View<int, 3> & v);
you may pass it an object of type ra::Owned<int, 3>.
Next: Rank extension, Previous: Containers and views, Up: Usage
To apply an operation to each element of an array, use the function for_each. The array is traversed in an order that is decided by the library.
ra::Small<double, 2, 3> a = {1, 2, 3, 4, 5, 6};
real s = 0.;
for_each([&s](auto && a) { s+=a; }, a);
⇒ s = 21.
To construct an array expression but stop short of executing it, use the function map. The expression will be executed implicitly when it is assigned to a view, printed out, etc.
using T = ra::Small<double, 2, 2>;
T a = {1, 2, 3, 4};
T b = {10, 20, 30, 40};
T c = map([](auto && a, auto && b) { return a+b; }, a, b); // (1)
⇒ c = { 11, 22, 33, 44 }
Expressions may take any number of arguments and be nested arbitrarily.
T d = 0;
for_each([](auto && a, auto && b, auto && d) { d = a+b; }, a, b, d); // same as (1)
for_each([](auto && ab, auto && d) { d = ab; },
map([](auto && a, auto && b) { return a+b; }, a, b), d); // same as (1)
The operator of an expression may return a reference and you may assign to an expression in that case. ra-ra will complain if the expression is somehow not assignable.
T d = 0;
map([](auto & d) -> decltype(auto) { return d; }, d) // just pass d along
= map([](auto && a, auto && b) { return a+b; }, a, b); // same as (1)
ra-ra defines many shortcuts for common array operations. You can of course just do:
T c = a+b; // same as (1)
Next: Subarray iteration, Previous: Array operations, Up: Usage
Rank extension is the mechanism that allows R+S to be defined even when R, S may have different ranks. The idea is an interpolation of the following basic cases.
Suppose first that R and S have the same rank. We require that the shapes be the same. Then the shape of R+S will be the same as the shape of either R or S and the elements of R+S will be
(R+S)(i₀ i₁ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ᵣ₋₁₎) + S(i₀ i₁ ... i₍ᵣ₋₁₎)
where r is the rank of R.
Now suppose that S has rank 0. The shape of R+S is the same as the shape of R and the elements of R+S will be
(R+S)(i₀ i₁ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ᵣ₋₁₎) + S().
The two rules above are supported by all primitive array languages, e.g. Matlab REF. But suppose that S has rank s, where 0<s<r. Looking at the expressions above, it seems natural to define R+S by
(R+S)(i₀ i₁ ... i₍ₛ₋₁₎ ... i₍ᵣ₋₁₎) = R(i₀ i₁ ... i₍ₛ₋₁₎ ... i₍ᵣ₋₁₎) + S(i₀ i₁ ... i₍ₛ₋₁₎).
That is, after we run out of indices in S, we simply repeat the elements. We have aligned the shapes so:
[n₀ n₁ ... n₍ₛ₋₁₎ ... n₍ᵣ₋₁₎] [n₀ n₁ ... n₍ₛ₋₁₎]
This rank extension rule is used by J REF and is known as prefix agreement. (The opposite rule of suffix agreement is used, for example, in Numpy REF.)
You can verify that the prefix agreement rule is distributive. Therefore it can be applied to nested expressions or to expressions with any number of arguments. This rule is applied systematically throughout ra-ra, even in assignment expressions. For example,
ra::Small<int, 3> x {3, 5, 9};
ra::Small<int, 3, 2> a = x; // assign x(i) to each a(i, j)
⇒ a = {{3, 3}, {5, 5}, {9, 9}}
ra::Small<int, 3> x(0.);
ra::Small<int, 3, 2> a = {1, 2, 3, 4, 5, 6};
x += a; // sum the rows of a
⇒ x = {3, 7, 11}
ra::Owned<double, 3> a({5, 3, 3}, ra::_0);
ra::Owned<double, 2> b({5}, 0.);
b += transpose<0, 1, 1>(a); // b(i) = ∑ⱼ a(i, j, j)
⇒ b = {0, 3, 6, 9, 12}
An obvious weakness of prefix agreement is that sometimes the axes you want to match are not the prefix axes. Obviously transposing the axes before extension is a workaround. For a more general form of rank extension, see The rank conjunction.
Other languages have a feature similar to rank extension called ‘broadcasting’. For example, in the way it’s implemented in Numpy, an array of shape [A B 1 D] will match an array of shape [A B C D]. The process of broadcasting consists in inserting so-called ‘singleton dimensions’ (axes with size one) to align the axes that one wishes to match. This is more general than prefix agreement, because any axes can be matched. You can think of prefix agreement as the case where all the singleton dimensions are added to the end of the shape.
Singleton broadcasting isn’t entirely general either, since it doesn’t handle transposition. Another drawback is that it muddles the distinction between a scalar and a vector of size 1. Sometimes, an axis of size 1 is no more than that, and if 2!=3 is a size error, it isn’t obvious why 1!=2 shouldn’t be.
Next: Slicing, Previous: Rank extension, Up: Usage
Next: Special objects, Previous: Subarray iteration, Up: Usage
Slicing is an array extension of the subscripting operation. However, tradition and convenience have given it a special status in most array languages, together with some peculiar semantics that ra-ra supports.
The form of the scripting operator A(i₀, i₁, ...) makes it plain that A is a function of rank(A) integer arguments. Therefore an array extension is immediately available through map. For example:
ra::Owned<double, 1> a = {1., 2., 3., 4.};
ra::Owned<int, 1> i = {1, 3};
map(a, i) = 77.;
⇒ a = {1., 77., 3, 77.}
Just as with any use of map, array arguments are subject to the prefix agreement rule.
ra::Owned<double, 2> a({2, 2}, {1., 2., 3., 4.});
ra::Owned<int, 1> i = {1, 0};
ra::Owned<double, 1> b = map(a, i, 0);
⇒ b = {3., 1.} // {a(1, 0), a(0, 0)}
ra::Owned<int, 1> j = {0, 1};
b = map(a, i, j);
⇒ b = {3., 2.} // {a(1, 0), a(0, 1)}
The latter is a form of sparse subscripting.
Most array operations (e.g. +) are defined through map in this way. For example, A+B+C is defined as map(+, A, B, C) (or the equivalent map(+, map(+, A, B), C)). Not so for the subscripting operation:
ra::Owned<double, 2> A({2, 2}, {1., 2., 3., 4.});
ra::Owned<int, 1> i = {1, 0};
ra::Owned<int, 1> j = {0, 1};
// {{A(i₀, j₀), A(i₀, j₁)}, {A(i₁, j₀), A(i₁, j₁)}}
ra::Owned<double, 2> b = A(i, j);
⇒ b = {{3., 4.}, {1., 2.}}
A(i, j, ...) is the outer product of the indices (i, j, ...) with operator A.
This operation sees much more use in practice than map(A, i, j ...). Besides, when the subscripts i, j, ... are scalars or linear ranges (integer sequences of the form (o, o+s, ..., o+s*(n-1)), the subscripting can be performed inmediately at constant cost and without needing to construct an expression object. This optimization is called beating REF. ra-ra is not smart enough to know when an arbitrary expression might be a linear range, so the following special objects are provided:
- Special object: iota count [start:0 [step:1]]
Create a linear range
start, start+step, ... start+step*(count-1).
This can used anywhere an array expression is expected.
ra::Owned<int, 1> a = ra::iota(4, 3 -2);
⇒ a = {3, 1, -1, -3}
Here, b and c are Views (see Containers and views).
ra::Owned<int, 1> a = {1, 2, 3, 4, 5, 6};
auto b = a(iota(3));
auto c = a(iota(3, 3));
⇒ a = {1, 2, 3}
⇒ a = {4, 5, 6}
- Special object: all
Create a linear range
0, 1, ... (nᵢ-1)when used as a subscript for thei-th argument of a subscripting expression.
This object cannot stand alone as an array expression. All the examples below result in View objects:
ra::Owned<int, 2> a({3, 2}, {1, 2, 3, 4, 5, 6});
auto b = a(ra::all, ra::all); // (1) a view of the whole of a
auto c = a(iota(3), iota(2)); // same as (1)
auto d = a(iota(3), ra::all); // same as (1)
auto e = a(ra:all, iota(2)); // same as (1)
auto f = a(0, ra::all); // first row of a
auto g = a(ra::all, 1); // second column of a
all is a special case (dots<1>) of the more general object dots.
- Special object: dots<n>
Equivalent to as many instances of
ra::allas indicated byn, which must not be negative. Each instance takes the place of one argument to the subscripting operation.
This object cannot stand alone as an array expression. All the examples below result in View objects:
auto h = a(ra::all, ra::all); // same as (1) auto i = a(ra::all, ra::dots<1>); // same as (1) auto j = a(ra::dots<2>); // same as (1) auto k = a(ra::dots<0>, ra::dots<2>); // same as (1) auto l = a(0, ra::dots<1>); // first row of a auto m = a(ra::dots<1>, 1); // second column of a
This is useful when writing rank-generic code, see examples/maxwell.C in the distribution for an example.
In addition to the special objects listed above, you can also omit any trailing ra::all subscripts. For example:
ra::Owned<int, 3> a({2, 2, 2}, {1, 2, 3, 4, 5, 6, 7, 8});
auto a0 = a(0); // same as a(0, ra::all, ra::all)
auto a10 = a(1, 0); // same as a(1, 0, ra::all)
⇒ a0 = {{1, 2}, {3, 4}}
⇒ a10 = {5, 6}
This supports the notion (see Rank polymorphism) that a 3-array is also an 2-array where the elements are 1-arrays themselves, or a 1-array where the elements are 2-arrays. This important property is directly related to the mechanism of rank extension (see Rank extension).
Next: The rank conjunction, Previous: Slicing, Up: Usage
Next: Compatibility, Previous: Special objects, Up: Usage
We have seen in Subarray iteration that it is possible to treat an r-array as an array of lower rank with subarrays as its elements. In the .iter<cell rank>() construction, this ‘exploding’ is performed (notionally) on the argument; the operation of the array expression is applied blindly to these cells, whatever they turn out to be.
for_each(sort, A.iter<1>()); // (in ra-ra) sort is a regular function, cell rank must be given for_each(sort, A.iter<0>()); // (in ra-ra) error, bad cell rank
The array language J has instead the concept of verb rank. Every function (or verb) has an associated cell rank for each of its arguments. Therefore .iter<cell rank>() is not needed.
for_each(sort_rows, A); // (not in ra-ra) will iterate over 1-cells of A, sort_rows knows
ra-ra doesn’t have ‘verb ranks’ yet. In practice one can think of ra-ra’s operations as having a rank of 0. However, ra-ra supports a limited form of J’s rank conjunction.
This is an operator that takes one verb (such operators are known as adverbs in J) and produces another verb with different ranks. These ranks are used for rank extension through prefix agreement, but then the original verb is used on the cells that result. The rank conjunction can be nested, and this allows repeated rank extension before the innermost operation is applied.
A standard example is ‘outer product’.
ra::Owned<int, 1> a = {1, 2, 3};
ra::Owned<int, 1> b = {40, 50};
ra::Owned<int, 2> axb = map(ra::wrank<0, 1>([](auto && a, auto && b) { return a*b; }),
a, b)
⇒ axb = {{40, 80, 120}, {50, 100, 150}}
It works like this. The verb ra::wrank<0, 1>([](auto && a, auto && b) { return a*b; }) has verb ranks (0, 1), so the 0-cells of a are paired with the 1-cells of b. In this case b has a single 1-cell. The frames and the cell shapes of each operand are:
a: 3 | b: | 2
Now the frames are rank-extended through prefix agreement.
a: 3 | b: 3 | 2
Now we need to perform the operation on each cell. The verb [](auto && a, auto && b) { return a*b; } has verb ranks (0, 0). This results in the 0-cells of a (which have shape ()) being rank-extended to the shape of the 1-cells of b (which is (2)).
a: 3 | 2 b: 3 | 2
This use of the rank conjunction is packaged in ra-ra as the from operator. It supports any number of arguments, not only two.
ra::Owned<int, 1> a = {1, 2, 3};
ra::Owned<int, 1> b = {40, 50};
ra::Owned<int, 2> axb = from([](auto && a, auto && b) { return a*b; }), a, b)
⇒ axb = {{40, 80, 120}, {50, 100, 150}}
Another example is matrix multiplication. For 2-array arguments C, A and B with shapes C: (m, n) A: (m, p) and B: (p, n), we want to perform the operation C(i, j) += A(i, k)*B(k, j). The axis alignment gives us the ranks we need to use.
C: m | | n A: m | p | B: | p | n
First we’ll align the first axes of C and A using the cell ranks (1, 1, 2). The cell shapes are:
C: m | n A: m | p B: | p n
Then we’ll use the ranks (1, 0, 1) on the cells:
C: m | | n A: m | p | B: | p | n
The final operation is a standard operation on arrays of scalars. In actual ra-ra syntax:
ra::Owned A({3, 2}, {1, 2, 3, 4, 5, 6});
ra::Owned B({2, 3}, {7, 8, 9, 10, 11, 12});
ra::Owned C({3, 3}, 0.);
for_each(ra::wrank<1, 1, 2>(ra::wrank<1, 0, 1>([](auto && c, auto && a, auto && b) { c += a*b; })), C, A, B);
⇒ C = {{27, 30, 33}, {61, 68, 75}, {95, 106, 117}}
Note that wrank cannot be used to transpose axes in general.
I hope that in the future something like C(i, j) += A(i, k)*B(k, j), where i, j, k are special objects, can be automatically translated to the requisite combination of wrank and perhaps also transpose. For the time being, you have to align or transpose the axes yourself.
Next: Extension, Previous: The rank conjunction, Up: Usage
Next: Internals, Previous: Compatibility, Up: Usage
For example:
ra::Owned<double, 1> x = map(cos, ra::Small<double, 1> {0.});
⇒ x = { 1. }
op can return a reference. A typical use is subscripting. For example (TODO better example, e.g. using STL types):
ra::Owned<int, 2> x = {{3, 3}, 0.};
map([](auto && i, auto && j) -> int & { return x(i, j); },
ra::Owned<int, 1> {0, 1, 1, 2}, ra::Owned<int, 1> {1, 0, 2, 1})
= 1;
⇒ x = {{0, 1, 0}, {1, 0, 1}, {0, 1, 0}}
Here the anonymous function can be replaced by simply x. Remember that unspecified return type defaults to (value) auto, so either a explicit type or decltype(auto) should be used to return a reference.
- Function: for_each op expr ...
Create array expression that applies op to expr ..., and traverse it.
op should normally return void. For example:
double s = 0.;
for_each([&s](auto && a) { s+=a; }, ra::Small<double, 1> {1., 2., 3})
⇒ s = 6.
- Function: pick select_expr expr ...
Create array expression that selects the first of expr ... if select_expr is 0, the second if select_expr is 1, and so on. The expressions that are not selected are not looked up.
This function cannot be defined using map, because map looks up each one of its argument expressions before calling op.
For example:
ra::Small<int, 3> s {2, 1, 0};
ra::Small<char const *, 3> z = pick(s, s*s, s+s, sqrt(s));
⇒ z = {1.41421, 2, 0}
- Function: where pred_expr true_expr false_expr
Create array expression that selects true_expr if pred_expr is
true, and false_expr if pred_expr isfalse. The expression that is not selected is not looked up.
For example:
ra::Owned<double, 1> s {1, -1, 3, 2};
s = where(s>=2, 2, s); // saturate s
⇒ s = {1, -1, 2, 2}
- Function: from op ... expr
Create outer product expression. This is defined as E(i00, i01 ..., i10, i11, ..., ...) = op(expr0(i01, i01 ...), expr1(i10, i11, ...) ...).
For example:
ra::Owned<double, 1> a {1, 2, 3};
ra::Owned<double, 1> b {10, 20, 30};
ra::Owned<double, 2> axb = from([](auto && a, auto && b) { return a*b; }, a, b)
⇒ axb = {{10, 20, 30}, {20, 40, 60}, {30, 60, 90}}
ra::Owned<int, 1> i {2, 1};
ra::Owned<int, 1> j {0, 1};
ra::Owned<double, 2> A({3, 2}, {1, 2, 3, 4, 5, 6});
ra::Owned<double, 2> Aij = from(A, i, j)
⇒ Aij = {{6, 5}, {4, 3}}
The last example is more or less how A(i, j) is actually implemented (see The rank conjunction).
- Function: wrank <input_rank ...> op
Wrap op using a rank conjunction (see The rank conjunction).
For example: TODO
⇒ x = 0
This operation does not work on arbitrary array expressions yet. TODO FILL
This is equivalent to view(ra::dots<k>, ra::iota(view.size(k), view.size(k)-1, -1)).
This operation does not work on arbitrary array expressions yet. TODO FILL
lo and hi are rank-1 expressions indicating the extent of the stencil on each dimension. Scalars are rank extended, that is, lo=0 is equivalent to lo=(0, 0, ..., 0) with length equal to the rank r of view. The stencil view has twice as many axes as view. The first r dimensions are the same as those of view except that they have their sizes reduced by lo+hi. The last r dimensions correspond to the stencil around each element of view; the center element is at s(i0, i1, ..., lo(0), lo(1), ...).
This operation does not work on arbitrary array expressions yet. TODO FILL
Next: The future, Previous: Reference, Up: Top
ra-ra applies its principles systematically and that can result in surprises.
behavior of operator= between views, construction of views.
a = b = c doesn’t work (expr = expr returns void).
rank extension in assignment of low rank = high rank.
- Item 0
- Item 1
- Item 2
Wishlist and acknowledged bugs.
You can still use pointer or std::initializer_list for shape by using the functions ptr or vector.