Numerical Arrays
Contents
Numerical Arrays
HIPP defines class templates for numerical arrays (i.e., vectors, matrices, or arrays of higher dimensionality). These types provide both expressive interface and efficient implementation, serving as the back-end of other numerical facilities.
Two favors of arrays are defined:
SArray: a static array whose rank and extents (i.e., shape) are determined at compile-time. Examples include a 3D position vector of a particle, a \(3\times 3\) inertial tensor of a rigid body, a histogram with a fixed number of 10 bins.DArray: a dynamic array whose rank is determined at compile-time, but the extent along each dimension are determined at run-time. Examples include a \(128\times 128\) single-channel images, a \(3\times 128 \times 128\) array that describes a 3 channel images, a \(512^3\) 3D spatial scalar field in numerical simulations.
The interface, underlying storage and performance are therefore different for those two types array.
In the following, we assume the declarations of namespace:
using namespace HIPP;
using namespace HIPP::NUMERICAL;
using namespace std;
which expose all the interface in the standard library and HIPP NUMERICAL module.
Static Arrays
Define SArray Instance
SArray class holds a static array. To get an instance, pass the scalar type of its element
and extents along dimensions as template parameters:
SArray<double, 3, 4> a34; // An array shaped [3][4].
SArray<double, 3> a3; // An vector shaped [3].
Because 1D arrays, i.e., vectors, are frequently used, we have a set of aliases:
SVec<t, N>class template: defines a vector with scalar typetandNelements. This is equivalent toSArray<t, N>.SVecXt<N>: a shortcut forSVec<t, N>. Heretcould bei(forint),f(forfloat),d(fordouble) andb(forbool).SVecNt: a shortcut forSVecXt<N>. HereNcould be1,2,3or4.
For example, the following three declarations all define the same vector type with 3 double-precision floating-point values:
SVec<double, 3> x3; // 3 double values.
SVecXd<3> y3;
SVec3d z3;
By default, all array elements are default-initialized (i.e., numerical values may have undefined values). To initialize the array elements at construction, simply use a initializer list with less of equal number of elements than the array:
SVec3d a3 = {1,2,3}, b3 = {1,2}; // b3[2] is uninitialized.
pout << "a3 = ", a3, ", b3 = ", b3, endl;
Output:
a3 = SArray{1,2,3}, b3 = SArray{1,2,2.31271e+35}
The following direct initializers are also acceptable:
double x[3] = {4,5,6};
SVec3d c3(1), // Constant array, {1,1,1}.
d3(x), // Copy from a pointer x.
e3(x, 2), // Copy only 2 elements from x.
f3(a3); // Copy from another SArray.
pout << "c3 = ", c3, ", d3 = ", d3, ", e3 = ", e3, ", f3 = ", f3, endl;
Output:
c3 = SArray{1,1,1}, d3 = SArray{4,5,6}, e3 = SArray{4,5,0}, f3 = SArray{1,2,3}
High-dimensional arrays can be initialized by treating them as 1D vector with row-major storage order:
SArray<int, 2, 3>
a23 = {1,2,3,4,5,6}, // Row-majored filling.
b23(1), // Constant array.
c23(&a23[0]), // Copy from a contiguous buffer.
d23(&a23[0], 3), // Copy only the first 3 elements.
e23(a23); // Copy from another SArray.
pout << "a23 = ", a23, ", b23 = ", b23, '\n',
"c23 = ", c23, ", d23 = ", d23, '\n',
"e23 = ", e23, endl;
Output:
a23 = SArray{1,2,3,4,5,6}, b23 = SArray{1,1,1,1,1,1}
c23 = SArray{1,2,3,4,5,6}, d23 = SArray{1,2,3,32576,-2005988080,21874}
e23 = SArray{1,2,3,4,5,6}
The value of the whole array can be changed later after construction:
a3 = c3; // Assign from another SArray.
b3 = 1.0; // Assign a scalar to all elements.
pout << "a3 = ", a3, ", b3 = ", b3, endl;
Output:
a3 = SArray{1,1,1}, b3 = SArray{1,1,1}
Element Access
SArray has tuple-like API defined for structured binding. For example, to bind
a list of variables to elements of a static vector, simply use:
SVec3i a3 = {1,2,3};
auto [x0,x1,x2] = a3; // Use structured binding.
pout << "x0,x1,x2 = ", x0, ",", x1, ",", x2, endl;
Output:
x0,x1,x2 = 1,2,3
Here, x0, x1, x2 are automatically declared as three int whose values
are initialized by the three elements of the vector.
Binding as reference is also allowed, by which the variables are declared as references to the vector elements. The update on the binding variables are reflected to the vector:
auto &[y0,y1,y2] = a3; // Bind as reference.
y0 += 1; y1 += 2; y2 += 3;
pout << "a3 = ", a3, endl;
Output:
a3 = SArray{2,4,6}
For higher dimensional array, the binding variables are declared as sub-arrays
along the first dimension of the parent array. For example, the following
code defines two vectors (i.e., SVec<int, 3>) whose values are initialized
by the two rows of the 2D array:
SArray<int, 2, 3> a23 = {1,2,3,4,5,6};
auto [z0, z1] = a23; // Bind to the 1st and 2nd rows.
pout << "z0 = ", z0, ", z1 = ", z1, endl;
Output:
z0 = SArray{1,2,3}, z1 = SArray{4,5,6}
Again, it is allowed to bind as reference:
auto &[w0, w1] = a23;
w0 += a3*10;
pout << "a23 = ", a23, endl;
Output:
a23 = SArray{21,42,63,4,5,6}
Like STL, SArray defines index-based and iterator-based element access
interface:
The method
data()returns a pointer to the first element of the array (always in a row-major order).size()returns the number of elements in the array.The method
raw()returns a reference to the raw array.begin()andend()returns the iterators.cbegin()andcend()are the constant counterparts.Operator
[]returns a reference to the element according to 1D index (row-major order). Methodat()is similar except that it throws an :class:ErrLogic`` exception on the out-of-range indexing.Operator
()accept N-Dimensional indices for a N-D array. This is the same as[]for 1D array, but more convenient for higher dimensionality.
For example, for a vector (1D array), its element access can be:
int *p3 = a3.data();
int (& a3_raw)[3] = a3.raw();
for(size_t i=0; i<a3.size(); ++i) {
p3[i] = i; // Visit the buffer by pointer.
a3_raw[i] += 1; // Visit by the raw array reference.
}
for(auto it = a3.begin(); it != a3.end(); ++it)
*it += 1; // Visit by the iterator.
a3[0] = -1; // Visit by operator [] indexing.
a3.at(1) = -2; // Visit by at().
pout << "a3 = ", a3, endl;
Output:
a3 = SArray{-1,-2,4}
For the 2D array:
int *p23 = a23.data();
int (&a23_raw)[2][3] = a23.raw();
p23[0] = -1; // (0,0) element.
a23_raw[0][1] = -2; // (0,1) element.
a23[2] = -3; // (0,2) element.
a23.at(3) = -4; // (1,0) element.
*(a23.begin() + 4) = -5; // (1,1) element.
a23(1,2) = -6; // (1,2) element.
pout << "a23 = ", a23, endl;
Output:
a23 = SArray{-1,-2,-3,-4,-5,-6}
To help interact with STL, 1D SArray has two methods, to_vector() and to_array,
which convert the instance to a STL vector and array, respectively:
vector<int> vec = a3.to_vector();
array<int, 3> arr = a3.to_array();
array<double, 3> arr_d = a3.to_array<double>(); // Change the target type.
pout << "vec = {", vec, "}, arr = {", arr, "}, arr<double> = {", arr_d, "}", endl;
Output:
vec = {-1,-2,4}, arr = {-1,-2,4}, arr<double> = {-1,-2,4}
Arithmetic Logic Operations
Element-wise arithmetic and logic operations are defined for SArray. Two operands
can be arrays of the same shapes, or an array with a scalar:
SVec3d a = {0,1,2}, b = {3,4,5};
pout << "a+1, a-1, a+b, a-b, a*b+10 => ",
a+1, ',', a-1, ',', a+b, ',', a-b, ',', a*b+10, endl;
Self-assignment versions of them are also available:
a += b;
b += 1;
pout << "a => ", a, ", b => ", b, endl;
Element-wise comparisions can also be made:
pout << "a==b, a!=b, a>0 => ",
a==b, ',', a!=b, ',', a>0, endl;
Output:
a+1, a-1, a+b, a-b, a*b+10 => SArray{1,2,3},SArray{-1,0,1},SArray{3,5,7},SArray{-3,-3,-3},SArray{10,14,20}
a => SArray{3,5,7}, b => SArray{4,5,6}
a==b, a!=b, a>0 => SArray{0,1,0},SArray{1,0,1},SArray{1,1,1}
Reduction
Reduction operations can be applied to one or more SArray instance to produce
scalar or other SArray instance.
norm(), squared_norm() and norm(p) returns the 2-norm (Frobenius norm), square of 2-norm, and p-norm
of an array:
SVec3d a = {0,1,2};
pout << "Norm = ", a.norm(),
", squared norm = ", a.squared_norm(),
", 1 norm = ", a.norm(1), endl;
normalize() perform in-place normalization on the instance. normalized() returns a normalized copy:
auto b = a.normalized(); // A normalized copy.
a.normalize(); // In-place normalization.
auto c = a.normalized(1); // Use 1 norm.
pout << "a = ", a, ", b = ", b, ", c = ", c, endl;
Output:
Norm = 2.23607, squared norm = 5, 1 norm = 3
a = SArray{0,0.447214,0.894427}, b = SArray{0,0.447214,0.894427}, c = SArray{0,0.333333,0.666667}
Minimum, maximum, mean, sum, product can be obtained by corresponding methods:
pout << "sum = ", a.sum(),
", prod = ", a.prod(),
", mean = ", a.mean(),
", min = ", a.min(),
", max = ", a.max(),
", min and max = ", a.minmax(),
", index of min = ", a.min_index(), endl;
Output:
sum = 1.34164, prod = 0, mean = 0.447214, min = 0, max = 0.894427, min and max = 0:0.894427, index of min = 0
any() tells whether or not the array has a non-false element, all() checks if all elements are non-false:
pout << "any true = ", a.any(),
", all true = ", a.all(),
", any positive = ", (a>0).any(),
", any equal = ", (a==b).any(), endl;
Output:
any true = 1, all true = 0, any positive = 1, any equal = 1
For 1D SArray, i.e., vector, dot() and cross() are defined for dot product and
cross product:
pout << "a dot b = ", a.dot(b),
", a cross b = ", a.cross(b), endl;
Output:
a dot b = 1, a cross b = SArray{0,0,0}
Here, because because a and b are normalized vectors and parallel to
each other, the dot product is 1 and the cross product is 0.
Apply Function to Array Elements
SArray defines both specific and general methods that map the instance element-wisely to a new one:
floor(),ceil(),trunc(): return a new integer instance whose elements are from the flooring, ceiling and truncation of the original elements, respectively. If applied to an integer array, they simply return a copy.abs(): returns a new instance with all elements to be the absolute values of original elements.mapped(unary_op): returns a new instance with each element to be the result of applying the callableunary_opto the original element.map(unary_op)is similar, but modify the instance in-place.visit(binary_op): call the callablebinary_op(i, x)for each elementxwith 1D indexi.
The following examples show how to use these methods:
SVec3d a = {-1.1, 1.1, 2.1};
pout << "floor = ", a.floor(),
", ceil = ", a.ceil(),
", trunc = ", a.trunc(),
", abs = ", a.abs(), endl;
Output:
floor = SArray{-2,1,2}, ceil = SArray{-1,2,3}, trunc = SArray{-1,1,2}, abs = SArray{1.1,1.1,2.1}
The methods mapped and map() enable more general transformation:
auto sin_a = a.mapped( [](auto x){ return sin(x); } ); // A mapped copy
a.map( [](double x) {return x+1.0; } ); // In-place change
pout << "sin_a = ", sin_a,
", a = ", a, endl;
Output:
sin_a = SArray{-0.891207,0.891207,0.863209}, a = SArray{-0.1,2.1,3.1}
The method visit() does not return a new instance nor modify the original instance. It simply
applies the user-passed function to the element and its index:
// Visit each element by index and value
a.visit( [](size_t i, auto x){ pout << '(', i, ')', x, '\n'; } );
Output:
(0)-0.1
(1)2.1
(2)3.1
View
HIPP provides “view” object which refers to a subset of elements in a SArray.
Two types of supported views are:
Boolean view: defined by a mask
SArraywhose extents are the same as the targetSArrayand whose elements arebool.Stride view: defined by a set of triplets in the form
(begin, end, step)at each dimension.
All the views are created by the view() method of SArray. For convenience,
a Boolean view can also be created by operator[] with a bool array passed as mask.
For example, the following code creates a mask array by comparison operator >
and then uses it to create a view:
SVec3d a = {-1, -2, -3};
auto pos_a = a[ a > 0 ]; // A view of all positive elements.
You may also explicitly create the mask and use it to filter elements:
SVec3b mask = {true, true, false};
auto first_2 = a[mask]; // A view of the first 2 elements.
View object has all arithmetic-logic operations defined as those for SArray.
Reduction operations sum(), prod(), mean(), all() and any(),
and visitors mapped(), map() and visit() are also defined.
When invoked, the corresponding operation is applied to only those selected
elements.
For example, use the above Boolean view:
first_2 += 2;
pout << a, endl;
The code changes only the first 2 elements of the parent SArray.
Output:
SArray{1,0,-3}
A stride view can be created by the view() method of any SArray, say,
a.
To ensure the overload exactly resolved, the first argument must be
a special object a.s_stride, and the following arguments must be
special slice object (created from) a.s_xxx. For example, to
select the range [0, 5) with step size 2, call:
SVec<double, 5> b = {1,2,3,4,5};
auto slice_b = b.view(b.s_stride, b.s_range(0,5,2)); // A view of the indexing range [0,5) with stride 2.
Then the view object can be used as:
slice_b += 10;
pout << b, endl;
Output:
SArray{11,2,13,4,15}
For high dimensional array, each dimension needs a slice object. For example, to select the first two rows of a 2D array, use:
SArray<int, 3, 3> c = {1,2,3,4,5,6,7,8,9};
auto slice_c = c.view(c.s_stride, c.s_head(2), c.s_all); // The first two rows
slice_c *= -1;
pout << c, endl;
Where c.s_head(n) select the first n elements at the corresponding dimension,
while c.s_all selects all.
Output:
SArray{-1,-2,-3,-4,-5,-6,7,8,9}