Chapter 8, part 9 of n: finite-difference iterators
Greetings.
This week, the third post in a series on the new finite-difference framework. It started here, in case you need to catch up.
Mandatory plug: there are still places availables for my Introduction to QuantLib Development course. For details, click on the link in the banner above.
Follow me on Twitter if you want to be notified of new posts, or add me to your Google+ circles, or subscribe via RSS: the buttons for that are in the footer. Also, make sure to check my Training page.
Operators, iterators and layouts
The interface of the FdmLinearOp
class is shown in the listing
below, and its main feature is an apply
method that returns the
result of applying the operator to a given array; you might remember a
similar applyTo
method defined by the old-style operators. (I’ll
gloss over the toMatrix
method, which returns some kind of
representation of the operator as a matrix and seems to be only used
for inspection.)
class FdmLinearOp {
public:
typedef Array array_type;
virtual ~FdmLinearOp() { }
virtual Disposable<array_type> apply(
const array_type& r) const = 0;
virtual Disposable<SparseMatrix> toMatrix() const = 0;
};
However, you probably noticed that the array_type
used in the
interface is defined as a typedef to Array
, which implements a
one-dimensional array. What of my claim that we can also cover
multi-dimensional operators, then? Well, the framework maps the
elements of a multi-dimensional array to the elements of a
one-dimensional one by flattening the former and turning it into the
latter.
As an example, let’s say that we have a three-dimensional array \( a_{i,j,k} \) with \( 0 \le i < M \), \( 0 \le j < N \), and \( 0 \le k < P \). We’ll map it into a one-dimensional sequence via a simple lexicographic ordering of the elements. We start from \( a_{0,0,0} \) as the first element, followed by \( a_{1,0,0} \) up to \( a_{M-1,0,0} \). Then we increase the second index and reset the first, obtaining \( a_{0,1,0} \) and following with \( a_{1,1,0} \) to \( a_{M-1,1,0} \). The next element is \( a_{0,2,0} \), and this repeats until we get to \( a_{M-1,N-1,0} \), at which point we follow with \( a_{0,0,1} \). You got the gist; we go through the whole thing again to \( a_{M-1,N-1,1} \), then \( a_{0,0,2} \) and so on until we finish with \( a_{M-1,N-1,P-1} \).
The FdmLinearOpIterator
and FdmLinearOpLayout
classes, whose
interfaces are shown below, cooperate to implement the logic above.
And yes, even though both classes are named after the FdmLinearOp
class, they actually work on the arrays to which operators are
applied.
class FdmLinearOpIterator {
public:
explicit FdmLinearOpIterator(const std::vector<Size>& dim);
explicit FdmLinearOpIterator(Size index = 0);
FdmLinearOpIterator(const std::vector<Size>& dim,
const std::vector<Size>& coordinates,
Size index)
void operator++();
bool operator!=(const FdmLinearOpIterator&);
};
class FdmLinearOpLayout {
public:
FdmLinearOpLayout(const std::vector<Size>& dim);
FdmLinearOpIterator begin() const;
FdmLinearOpIterator end() const;
Size index(const std::vector<Size>& coordinates) const;
Size neighbourhood(const FdmLinearOpIterator& iterator,
Size i, Integer offset) const;
Size neighbourhood(const FdmLinearOpIterator& iterator,
Size i1, Integer offset1,
Size i2, Integer offset2) const;
};
Instances of the FdmLinearOpIterator
keep track of a tuple of
indices into the several dimensions of an array (like the \( (i,j,k) \)
tuples above) and of the corresponding single index into the flattened
array; to do that, they must also store the sizes of the array along
each dimension (the ones that I called \( M \), \( N \) and \( P \)
in the example).
Depending on the constructor you’ll call, you can pass explicitly all
the state or just part of it. The constructor taking only a vector of
dimensions stores the sizes of the array and implicitly initializes
all the indices to 0; that is, it returns the iterator pointing to the
beginning of the array. The constructor taking a single index builds
a kind of sentinel value; it cannot be incremented, because it doesn’t
know the size of the array, but when it is passed the total number of
elements of the array it returns an end
iterator to which
others can be compared to end a loop, as in:
std::vector<Size> dims(...); // { M, N, P }
FdmLinearOpIterator begin(dims);
FdmLinearOpIterator end(M*N*P);
for (FdmLinearOpIterator i=begin; i!=end; ++i)
... // do something
Finally, the third constructor takes the dimensions and the whole set of indices, and thus allows you to create an iterator pointing to any element of the array. For performance reasons (I guess) it leaves to the caller to check that the parameters are consistent, that is, that the set of indices into the array actually corresponds to the single index into the flattened array.
A FdmLinearOpIterator
instance is not really an iterator in the C++
sense, since it doesn’t implement the whole required interface; in
particular, there’s no dereference operator returning the value it
points to. What it does implement is operator++
, which increments
the iterator, and operator!=
, which compares two iterators; the two
methods are shown in the next listing. (Providing
only operator!=
and not operator==
is enough to support the i !=
end
idiom used in for loops, but might be a tad extreme.)
void FdmLinearOpIterator::operator++() {
++index_;
for (Size i=0; i < dim_.size(); ++i) {
if (++coordinates_[i] == dim_[i])
coordinates_[i] = 0;
else
break;
}
}
bool FdmLinearOpIterator::operator!=(
const FdmLinearOpIterator& iterator) {
return index_ != iterator.index_;
}
To round up their interface, iterators also declare the coordinates
inspector, which returns the set of indices into the array, and
index
, returning the corresponding single index.
The increment operator literally implements the iterating strategy I
described in the previous page. It starts by incrementing the first
index; if doing so causes it to reach the size of the first dimension,
it resets the index to 0 and repeats both the increment and check on
the next dimension; otherwise, it breaks out of the loop. Thus, for
instance, calling operator++
with the indices set to \( (M-1,N-1,2)
\) will increment the first and get to \( (M,N-1,2) \), then reset
it and increment the second to \( (0,N,2) \), then reset the second
too and increment the third to \( (0,0,3) \), then finally exit.
Calling it with the indices set to \( (2,0,0) \), instead, will just
increment the first to \( (3,0,0) \) and exit. In any case, the
single index into the flattened array is also incremented so that it
stays in sync with the other indices.
Finally, operator!=
just compares two iterators by comparing
their two indices into the flattened array (also known as the indices
for which I should probably get a shorter name at this point). This
allows iterators to work as sentinels even if they don’t have a full
set of indices and can’t be incremented.
The FdmLinearOpLayout
class (whose interface we have seen in a
previous listing, and whose implementation is shown in the next
listing) provides some utilities for working with iterators based on a
given array layout, or set of dimensions.
FdmLinearOpLayout::FdmLinearOpLayout(
const std::vector<Size>& dim)
: dim_(dim), spacing_(dim.size()) {
spacing_[0] = 1;
std::partial_sum(dim.begin(), dim.end()-1,
spacing_.begin()+1, std::multiplies<Size>());
size_ = spacing_.back()*dim.back();
}
FdmLinearOpIterator FdmLinearOpLayout::begin() const {
return FdmLinearOpIterator(dim_);
}
FdmLinearOpIterator FdmLinearOpLayout::end() const {
return FdmLinearOpIterator(size_);
}
Size FdmLinearOpLayout::index(
const std::vector<Size>& coordinates) const {
return std::inner_product(coordinates.begin(),
coordinates.end(),
spacing_.begin(), Size(0));
}
Most of its calculations are based on what it calls spacing between the elements of the array along the various directions. I’ll use the same example as before, and consider a three-dimensional array with sizes \( (M,N,P) \) along the three dimensions. If we take two elements \( (i,j,k) \) and \( (i+1,j,k) \), they occupy consecutive places on the one-dimensional flattened array (remember that we increment the first index first while enumerating the elements) so we say that the spacing between them is \( 1 \). If we take two elements \( (i,j,k) \) and \( (i,j+1,k) \), their positions in the flattened array are separated by a whole turn of the first index: enumerating the elements between them will cause you to go through \( (i+1,j,k) \), \( (i+2,j,k) \) and so on up to \( (M-1,j,k) \), then to \( (0,j+1,k) \), and then again through \( (1,j+1,k) \), \( (2,j+1,k) \) up to \( (i,j+1,k) \). A bit of thinking will convince you that the spacing between the two elements is \( M \), and a similar argument will find a spacing of \( M \times N \) between elements \( (i,j,k) \) and \( (i,j,k+1) \).
The constructor of a FdmLinearOpLayout
instance takes a general set
of sizes \( (n_1, \ldots, n_N) \) and calculates the spacings \(
(1, n_1, n_1 \times n_2, \ldots, \prod_1^{N-1} n_i) \) as well as the
total size \( \prod_1^N n_i \); the calculations is written in STL
speak, but it shouldn’t be difficult to map it to the formulas.
Given the sizes and the spacings, the class can implement a number of
methods that let you (and library code) work with iterators. The
begin
and end
methods return the corresponding
iterators, so that once you have a layout you no longer need to care
about sentinels and iterator constructors: loops
can be written as:
std::vector<Size> dims(...);
FdmLinearOpLayout l(dims);
for (FdmLinearOpIterator i=l.begin(); i!=l.end(); ++i)
... // do something
The index
method provides a conversion between a set of indices
and the single index into the flattened array. Other methods, such as
neighborhood
, work at a slightly higher level; they return the
index (or the iterator) of the element you get by starting from a
given element, whose position is described by an iterator, and moving
a given number of steps in a given dimension. An overload that moves in
two directions is also available for convenience, even though it’s
equivalent to two consecutive calls to the first version.
At this point, we can finally make sense of the interface of the
FdmMesher
class shown in the previous post. The methods dplus
, dminus
and location
all take as first argument an iterator to specify a
point in the mesh; the second argument specifies the direction along
which we’re moving. (The points in the mesh are in one-to-one
correspondence with the elements of the array containing the values of
the instrument being priced on the mesh, thus an iterator can denote
either one.) Finally, the layout
method returns an instance of
FdmLinearOpLayout
that can be used to create iterators.
In next post: a few actual operators.