Tutorial

Environment

Dependencies

The versions below indicate the libraries !ndarray is currently developed and tested with. New minor versions should generally work as well, but have not been tested. ndarray makes extensive use of C++ template metaprogramming, and may not work with older or non-standard-compliant compilers.

• Core ndarray library (ndarray.h): boost 1.38-1.42
• Python conversion module (ndarray/python.h): python 2.6, numpy 1.2
• Boost.Python module (ndarray/python/boost/): boost.python 1.38-1.42, python 2.6, numpy 1.2
• Eigen interface (ndarray/eigen.h): Eigen 2.0
• Fast Fourier transforms (ndarray/fft.h): FFTW 3.2

Installation

ndarray is a header-only library; after downloading and unpacking the source, you can start using it immediately simply by adding it to your compiler's include path.

For tests, we use the SCons build system, but SCons is not necessary to make use of the library.

Creating New Arrays

Array objects have two normal constructors intended for public use, the default constructor and a converting copy constructor. The default constructor creates an "empty" array, with a null data pointer and zero shape. The copy constructor creates a shared view into an existing array.

New Memory

To create a new non-trivial array, one can use the allocate function, which returns a temporary object implicitly convertible to Array:

* Array<double,3,3> a = allocate(makeVector(3,4,5));
* // equivalent to
* // >>> a = numpy.empty((3,4,5),dtype=float)
* // in Python
* 

The makeVector function here is a variable-argument-length constructor for the Vector object, a fixed-size array class whose int variant is used to specify shape and strides for Array. The appropriate Vector template for a particular Array template is available as the Index typedef within the Array class.

The allocate function can also take an STL allocator as a template argument and/or regular argument:

* Array<double,3,3> a = allocate< std::allocator<void> >(makeVector(3,4,5));
* 

* Array<double,3,3> a = allocate(makeVector(3,4,5), std::allocator<void>());
* 

Note that the type of the allocator does not have to match the type of the array; the allocator's "rebind" member will be used to create the correct allocator when the array is constructed. Furthermore, unlike standard library containers, Array is not templated on its allocator type; after construction, it is impossible to determine how an Array's memory was allocated. An Array constructed by allocate is not generally initialized to any value (do not assume it contains zeros).

External Memory

Arrays can also be constructed that point to external data:

* #include <cstdlib>
* Array<double,1,1>::Owner owner(std::malloc(sizeof(double)*5), std::free);
* Array<double,1,1>::Index shape = makeVector(5);
* Array<double,1,1>::Index strides = makeVector(1);
* Array<double,1,1> a = external(owner.get(), shape, strides, owner);
* 

The 'strides' vector here specifies the space between elements in each dimension; the dot product of the strides vector with an index vector should give the offset of the element with that index from the first element of the array. The 'Owner' type here is a typedef to a boost::shared_ptr, which can take an arbitrary functor as a custom deleter (here, std::free). By defining an appropriate deleter, an array can manage virtually any kind of memory. However, it is also possible to create an array with no reference counting by passing an empty owner (or passing none at all):

* #include <cstdlib>
* double data[] = { 5.3, 1.2, 6.3, 2.8, 7.0 };
* Array<double,1,1>::Index shape = makeVector(5);
* Array<double,1,1>::Index strides = makeVector(1);
* Array<double,1,1> a = external(data, shape, strides);
* 

In this case, the user is responsible for ensuring that the data pointer provided to the array remains valid during the array's lifetime, and is eventually deallocated later.

Assignment

Direct Array assignment is shallow:

* Array<double,1,1> a = allocate(makeVector(5));
* Array<double,1,1> b;
* b = a; // the two arrays now share data.
* 

To actually set the elements of an array, we can use Array::deep():

* Array<double,11,> b = allocate(a.getShape());
* b.deep() = a; // 'b' is now a deep copy of 'a'.
* 

Scalar assignment and augmented assignment operations are also supported:

* b.deep() = 5.2;
* b.deep() += a;
* 

The deep() method returns a proxy ArrayRef object, which behaves just like an Array aside from its assignment operators, and is implicitly convertible to Array.

Indexing and Iteration

A multidimensional Array behaves like a container of Arrays with lower dimensions, while a one-dimensional Array behaves like a container of elements:

* int data = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
* Array<int,2,2> a = external(data, makeVector(4,3));
* // make 'b' a view into the second row of 'a' (Reference is a typedef to Array<int,1,1>):
* Array<int,2,2>::Reference br = a[1]; // br is an ArrayRef<int,1,1>.
* Array<int,2,2>::Value b = a[1];      // b is an Array<int,1,1>
* Array<int,2,2>::Element b0 = b[0];   // b0 == 6; Element is a typedef to int.
* Array<int,2,2>::Reference::Value b1 = b[1];     // b1 == 7; Reference::Value is also int.
* Array<int,2,2>::Reference::Reference b2 = b[2]; // b2 == 8; Reference::Reference is int&.
* 

Indexing operations return ArrayRef objects, not Arrays. This allows them to be assigned to without manually calling deep():

* a[1] -= 3; // subtract three from the entire second row.
* 

For one dimensional arrays, the "Reference" typedef is equivalent to "Element &", while the "Value" typedef is equivalent to "Element". For multidimensional arrays, "Reference" is a lower-dimensional ArrayRef, while "Value" is a lower-dimensional Array.

Array is designed to have lightweight nested iterators, with types provided by the Array::Iterator typedef. For contiguous one-dimensional arrays (Array<T,1,1>), this is a typedef to a simple pointer. The typical pattern to iterate over a 3-dimensional array looks like the following:

* Array<double,3,3> a = allocate(makeVector(5,6,8));
* for (Array<double,3,3>::Iterator i = a.begin(); i != a.end(); ++i) {
*     for (Array<double,3,3>::Reference::Iterator j = i->begin(); j != i->end(); ++j) {
*         for (Array<double,3,3>::Reference::Reference::Iterator k = j->begin(); k != j->end(); ++k) {
*              // *k == a[i - a.begin()][j - i->begin()][k - j->begin()];
*         }
*     }
* }
* 

As expected, the iterators of multidimensional arrays dereference to lower-dimensional arrays, and the iterators of one-dimensional arrays dereference to elements. With some compilers, it may be advantageous to move the calls to end() outside their loops.

Just like direct indexing, multidimensional array iterators dereference to ArrayRef, not Array.

STL-compliant typedefs "iterator", "const_iterator", "reference", "const_reference", and "value" are also provided, though the const variants are not actually const (because Array provides smart-pointer const-correctness rather than container const-correctness).

Single elements can be extracted from multidimensional arrays by indexing with ndarray::Vector:

* a[makeVector(3,2,1)] == a[3][2][1];
* 

Views

General views into a an Array are created by passing a ViewDef object to the [] operators of Array, returning a new Array that shares data and owns a reference to the original.

ViewDef involves a lot of template metaprogramming, so the actual template class is an internal detail, and ViewDefs are constructed by calls to the view() function function followed by chained calls to the function call operator, resulting in a syntax that looks like this:

* Array<double,5> a = allocate(makeVector(3,5,2,6,4));
* b = a[view(1)(0,3)()(4)];
* 

which is equivalent to the Python code:

* a = numpy.empty((3,5,2,6,4),dtype=float)
* b = a[1,0:3,:,4]
* 

The value passed to each call specifies how to extract values from that dimension:

• A single integer selects a single subarray from that dimension, reducing the overall dimensionality of the view relative to the parent array by one.
• An empty call selects the entire dimension.
• A pair of integers selects a contiguous subset of subarrays from that dimension.
• A triplet of integers selects a noncontiguous subset of subarrays from that dimension.

Any dimensions which are not specified because the length of the ViewDef expression is smaller than the dimensionality of the parent array will be considered full-dimension selections:

* a[view(3)] == a[view(3)()()()()];
* 

Arithmetic Operators and Comparison

Arrays provide element-wise arithmetic operations that make use of expression templates:

* Array<double,2,2> a = allocate(makeVector(3,4));
* // initialize the elements of 'a'
* Array<double,2,2> b = allocate(a.getShape());
* b = a * 3 + 2; // this expression is lazy, and never allocates a temporary array
* 

We can simplify the previous example by initializing 'b' with the copy() function, which is simply a shortcut for allocate and assign:

* Array<double,2,2> a = allocate(makeVector(3,4));
* // initialize the elements of 'a'
* Array<double,2,2> b = copy(a * 3 + 2);
* 

As a rule, ndarray never allocates memory for a new unless you explicitly tell it to with the allocate() or copy() functions.

Element-wise comparisons are also supported, but not via overloaded operators:

* Array<double,2,2> a = allocate(makeVector(3,4));
* // initialize the elements of 'a'
* Array<double,2,2> b = allocate(makeVector(3,4));
* // initialize the elements of 'b'
* Array<bool,2,2> c = copy(equal(a, b));
* 

The element-wise comparison functions (equal, not_equal, less, less_equal, greater, greater_equal) and logical operators (logical_and, logical_or, logical_not) are also lazy, and are most useful when used in conjunction with the reduction functions all() and any():

* Array<double,2,2> a = allocate(makeVector(3,4));
* // initialize the elements of 'a'
* Array<double,2,2> b = allocate(makeVector(3,4));
* // initialize the elements of 'b'
* bool v1 = all(logical_or(greater(a, b), greater(a, 3.0)));
* bool v2 = any(less(b, a));
* 

Array does overload the equality and inequality operators, but these compare for "shallow" equality, not element-wise equality:

* Array<double,2,2> a = allocate(makeVector(3,4));
* Array<double,2,2> b = copy(a);
* bool v1 = (a == b); // false, because 'a' and 'b' do not share data.
* Array<double const,2,1> c(a);
* bool v2 = (a == c); // true, because 'a' and 'c' have the same data, shape, and strides.
*