ra-ra/examples

// (c) Daniel Llorens - 2016

// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 3 of the License, or (at your option) any
// later version.

// @file readme.C
// @brief Examples used in top-level README.md

// TODO Generate README.md and/or these examples.

#include "ra/operators.H"
#include "ra/io.H"
#include "ra/test.H"
#include

using std::cout; using std::endl;

int main()
{
TestRecorder tr(std::cout);
tr.section("dynamic/static shape");
// Dynamic or static array rank. Dynamic or static array shape (all dimensions or none).
{
ra::Big A({2, 3}, 'a'); // dynamic rank = 2, dynamic shape = {2, 3}
ra::Big B({2, 3}, 'b'); // static rank = 2, dynamic shape = {2, 3}
ra::Small C('c'); // static rank = 2, static shape = {2, 3}
cout << "A: " << A << "\n\n";
cout << "B: " << B << "\n\n";
cout << "C: " << C << "\n\n";
}
tr.section("storage");
// Memory-owning types and views. You can make array views over any piece of memory.
{
// memory-owning types
ra::Big A({2, 3}, 'a'); // storage is std::vector inside A
ra::Unique B({2, 3}, 'b'); // storage is owned by std::unique_ptr inside B
ra::Small C('c'); // storage is owned by C itself, on the stack

cout << "A: " << A << "\n\n";
cout << "B: " << B << "\n\n";
cout << "C: " << C << "\n\n";

// view types
char cs[] = { 'a', 'b', 'c', 'd', 'e', 'f' };
ra::View D1({2, 3}, cs); // dynamic sizes and strides, C order
ra::View D2({{2, 1}, {3, 2}}, cs); // dynamic sizes and strides, Fortran order.
ra::SmallView, mp::int_list<3, 1>> D3(cs); // static sizes & strides, C order.
ra::SmallView, mp::int_list<1, 2>> D4(cs); // static sizes & strides, Fortran order.

cout << "D1: " << D1 << "\n\n";
cout << "D2: " << D2 << "\n\n";
cout << "D3: " << D3 << "\n\n";
cout << "D4: " << D4 << "\n\n";
}
tr.section("shape agreement");
// Shape agreement rules and rank extension (broadcasting) for rank-0 operations of any arity
// and operands of any rank, any of which can a reference (so you can write on them). These
// rules are taken from the array language, J.
// (See examples/agreement.C for more examples.)
{
ra::Big A({2, 3}, { 1, 2, 3, 1, 2, 3 });
ra::Big B({-1, +1});

ra::Big C({2, 3}, 99.);
C = A * B; // C(i, j) = A(i, j) * C(i)
cout << "C: " << C << "\n\n";

ra::Big D({2}, 0.);
D += A * B; // D(i) += A(i, j) * C(i)
cout << "D: " << D << "\n\n";
}
tr.section("rank iterators");
// Iterators over cells of arbitrary rank.
{
ra::TensorIndex<0> i;
ra::TensorIndex<1> j;
ra::TensorIndex<2> k;
ra::Big A({2, 3, 4}, i+j+k);
ra::Big B({2, 3}, 0);
cout << "A: " << A << "\n\n";

// store the sum of A(i, j, ...) in B(i, j). All these are equivalent.
B = 0; B += A; // default agreement matches prefixes
for_each([](auto && b, auto && a) { b = ra::sum(a); }, B, A); // default agreement matches prefixes
for_each([](auto && b, auto && a) { b = ra::sum(a); }, B, A.iter<1>()); // give cell rank
for_each([](auto && b, auto && a) { b = ra::sum(a); }, B, A.iter<-2>()); // give frame rank
cout << "B: " << B << "\n\n";

// store the sum of A(i, ...) in B(i, j). The op is re-executed for each j, so don't do it this way.
for_each([](auto && b, auto && a) { b = ra::sum(a); }, B, A.iter<2>()); // give cell rank
cout << "B: " << B << "\n\n";
}
// A rank conjunction (only for static rank and somewhat fragile).
tr.section("rank conjuction");
{
// This is a translation of J: A = (i.3) -"(0 1) i.4, that is: A(i, j) = b(i)-c(j).
ra::Big A = map(ra::wrank<0, 1>(std::minus()), ra::iota(3), ra::iota(4));
cout << "A: " << A << "\n\n";
}
// A proper selection operator with 'beating' of range or scalar subscripts.
// See examples/slicing.C for more examples.
tr.section("selector");
{
// TODO do implicit reshape in constructors?? so I can accept any 1-array and not only an initializer_list.
ra::Big A({2, 2, 2}, {'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h'});
cout << "A: " << A << "\n\n";

// these are all equivalent to e.g. A(:, 0, :) in Octave.
cout << "A1: " << A(ra::all, 0) << "\n\n";
cout << "A2: " << A(ra::all, 0, ra::all) << "\n\n";
cout << "A3: " << A(ra::all, 0, ra::dots<1>) << "\n\n";

// an inverted range.
cout << "A4: " << A(ra::iota(2, 1, -1)) << "\n\n";

// indices can be arrays of any rank.
ra::Big I({2, 2}, {0, 3, 1, 2});
ra::Big B({4}, {'a', 'b', 'c', 'd'});
cout << "B(I): " << B(I) << "\n\n";

// multiple indexing performs an implicit outer product. this results in a rank
// 4 array X = A(J, 1, J) -> X(i, j, k, l) = A(J(i, j), 1, J(k, l))
ra::Big J({2, 2}, {1, 0, 0, 1});
cout << "A(J, 1, J): " << A(J, 1, J) << "\n\n";

// explicit indices do not result in a View view (= pointer + strides), but the
// resulting expression can still be written on.
B(I) = ra::Big({2, 2}, {'x', 'y', 'z', 'w'});
cout << "B: " << B << endl;
}
// A TensorIndex object as in Blitz++ (with some differences).
tr.section("tensorindex");
{
// as shown above.
}
tr.section("STL compat");
{
ra::Big A = {'x', 'z', 'y'};
std::sort(A.begin(), A.end());
cout << "A: " << A << "\n\n";

ra::Big B({2, 2}, {1, 2, 3, 4});
B += std::vector({10, 20});
cout << "B: " << B << "\n\n";
}
// example from the manual [ma100].
{
ra::Small s {2, 1, 0};
ra::Small z = pick(s, s*s, s+s, sqrt(s));
cout << "z: " << z << endl;
}
// example from the manual [ma101].
{
ra::Big A({2, 5}, "helloworld");
std::cout << format_array(transpose<1, 0>(A), false, "|") << std::endl;
}
{
ra::Big A = {"hello", "array", "world"};
std::cout << format_array(A, false, "|") << std::endl;
}
// example from the manual [ma102].
{
// ra::Big A({3}, "hello"); // ERROR b/c of pointer constructor
ra::Big A({3}, ra::scalar("hello"));
std::cout << format_array(A, false, "|") << std::endl;
}
// example from the manual [ma103].
{
ra::Big A({3, 2}, {1, 2, 3, 4, 5, 6});
ra::Big B({2, 3}, {7, 8, 9, 10, 11, 12});
ra::Big 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);
/* 3 3
27 30 33
61 68 75
95 106 117 */
cout << C << endl;
}
// example from the manual [ma104]
{
ra::Big c({3, 2, 2}, ra::_0 - ra::_1 - 2*ra::_2);
cout << c << endl;
cout << map([](auto && a) { return sum(diag(a)); }, iter<-1>(c)) << endl;
}
// example from the manual [ma105]
{
ra::Big a({2, 3}, {1, 2, 3, 4, 5, 6});
ra::Big b({3}, {10, 20, 30});
ra::Big c({2, 3}, 0);
iter<1>(c) = iter<1>(a) * iter<1>(b); // multiply each item of a by b
cout << c << endl;
}
// example from the manual [ma109]. This is a rare case where I need explicit ply.
{
ra::Big o = {};
ra::Big e = {};
ra::Big n = {1, 2, 7, 9, 12};
ply(where(odd(n), map([&o](auto && x) { o.push_back(x); }, n), map([&e](auto && x) { e.push_back(x); }, n)));
cout << "o: " << format_array(o, false) << ", e: " << format_array(e, false) << endl;
}
return 0;
}