Usage of Boost and STL in Scientific Programming

Tips and Tricks

Code snippet below is a function that is used in one of our module (programming environment) called DamBreak++ Simulator (link on the web page).

// initialize the field with some initial values
// for our test, we use the dam-break problem
void DamBreakInit( TestField& aField)
{
// namespace
using namespace std;
using namespace boost::lambda;

// retrieve the spatial discretization
TestField::DiscrInfo w_Disrc = aField.getDiscrInfo();

// initializer (unary functor) which set initial values
DamBreakStepFunc w_testForeachptr(w_Disrc);

// these 2 below work perfectly with standard algorithm
// initialization of nodal values (field)
for_each( aField.values().begin(), aField.values().end(), w_testForeachptr);

// another approach by using the bind mechanism
for_each( aField.begin(), aField.end(), bind(Evaluation_A_fonction_H, _1));
}

The method argument “Field” is a type that hold global nodal variables stored in a pointer container (boost pointer container).

class TestField
{
public:
typedef std::shared_ptr<Nvalue> vaptr;
typedef boost::ptr_vector<Nvalue>::value_type value_type;
typedef boost::ptr_vector<Nvalue>::iterator iterator;
typedef boost::ptr_vector<Nvalue>::const_iterator const_iterator;
public:
// …
iterator begin() { return m_NodeValues.begin();}
iterator end() { return m_NodeValues.end(); }
const_iterator begin() const { return m_NodeValues.begin();}
const_iterator end() const { return m_NodeValues.end(); }

…..

// return write/read access (deprecated)
boost::ptr_vector<Nvalue>& values() { return m_NodeValues;}
const boost::ptr_vector<Nvalue>& values() const { return m_NodeValues;}

….

These nodal variables (“NValue”) represent state and other variables at grid node and are updated (modified) as the simulation step through time. We need a reference semantic and pointer container offers this functionality (more specifically the pointer vector).
Reminder: main advantage of using boost pointer container instead of standard vector of pointer, it let you treat pointer container exactly the same way as a container of objects.

First call to “for_each” algorithm

The first 2 lines of the “DamBreakInit” function retrieve the discretization parameters for the DamBreak solution and create a function (math).

// retrieve the spatial discretization
TestField::DiscrInfo w_Disrc = aField.getDiscrInfo();

// initializer (unary functor) which set initial values
DamBreakStepFunc w_testForeachptr(w_Disrc);

DamBreakStepFunc” is a function that initialize the Dam Break problem which is represented by a mathematical step function. Actually this function is called a “functor”. Functor object For those who are not familiar with functor, a “functor” (or functional, functionoid, or simply function object) is an object that masquerades as a function, it means that you can call it like ‘f(x)’. In scientific computing is just more natural to use such objects which are close to real world type (math function). Below the header of the function that inherits from std::unary_function.

class DamBreakStepFunc :
public std::unary_function< TestBoost::Nvalue, void>
{
public:
// default ctor
DamBreakStepFunc( TestField::DiscrInfo aDiscr)
: m_Discr( aDiscr) // at the beginning
{}

// to be used with the for_each algorithm
result_type operator() (argument_type& aNval) const
{do some stuff}

It is important to note that we could have just declared a “functor” without inheriting from std::unary_function, but by doing this, we wouldn’t have inherited some structures from base type (unary_function) such as the “result_type” and “argument_type” that are needed by many of the STL algorithm. It is a good practice to inherit from the unary function. Last thing, we need to override the operator “()” since by definition of a “functor” has the following signature (expected by STL algorithm too). Here “argument_type” is a reference to an object of the collection and “result type” is the return type of the function, in our case it’s void, no return.

Then we use the “for_each” from stl algorithm to set variables initial values.

// initialization of nodal values (field)
for_each( aField.values().begin(), aField.values().end(), w_testForeachptr);

“for_each” algorithm shall be used whenever is possible to avoid using for-loop which are prone to errors (index initialization and so on). Moreover, “for_each” state it clear; we are going through each element of the collection, the range [begin,end], by applying the unary function (function with only one argument) also called “functor” to each object. First 2 arguments of the “for_each” algorithm define the range on which we apply the function and is expressed by “aField.values().begin()/end()” pair. It returns an iterator from pointer container which is a reference to the object (pointer container iterators iterate over pointed-to objects, not pointers, see my blog on an introduction on boost pointer container “Basics of Pointer Container: Part I”). Moreover same syntax as with standard container like vector (that’s make it so beautiful). The last argument stand for the function objects that represent a step function.

 Second call to “for_each” algorithm call (C++11 feature)

// using bind mechanism
for_each( aField.begin(), aField.end(), bind(Evaluation_A_fonction_H, _1));

We use again the “for_each” algorithm but this time we make use of the new feature of C++11 the bind mechanism.  “Evaluation_A_fonction_H” in the bind argument is a free function which take as argument the NValue (nodal value).

void Evaluation_A_fonction_H( Nvalue& aNodeVal)
{
aNodeVal[0]=aNodeVal[2]; for example
}

Bind mechanism was introduced by boost library and then moved in C++1 whose main purpose is creating function objects on-the-fly. When using algorithm from the stl, we often need to supply function object  such as for_each (above). But sometimes we don’t want to write a function object, we just want to use some free function or member function and pass it to the algorithm (customize). STL provide adaptor such as bind1st and bind2nd, but it is limited on what you can do. That’s where bind comes to play. Conceptually, bind is a generalization of the existing standard library function bind1st and bind2nd, with additional functionality that allows for more sophisticated functional composition.

bind mechanism The code above illustrates such a class, wrapping a function pointer up in an object by using the bind adaptor. The adaptor above is an example of the Adapter Pattern, whose intent is described in Design Patterns as

Convert the interface of a class into another interface clients expect. Adapter lets classes work together that couldn’t otherwise because of incompatible interfaces”.

 Conclusion

Using Boost and STL (C++11) bring the mathematical programming to new level. First, write less code and easier to understand end eventually maintain. Scientific programming involves loop on a set of data and generally written with for-loop which, as we have already mentioned, error prone and difficult to understand. When working  in a “Test Driven Development” each component is tested, In our case, this mean the functorDamBreakStepFunc” has been already tested and if there an error occur it has been already flagged (suppose that you wrote good unit test).

clients and partners

Autolog