Template Meta Programming Techniques

[Jon Rea]

I work with molecular simulations. In a normal infinate space
simulation, a particle can occupy any cartesian position in space.

If i want the distance between two particles in a system I can use:

Maths::dvector va, vb;
// Code here: Do something to initialise variables and move them around...
double distance = dist(va,vb);

I am currently looking at 'Periodic boundary simulations'. By this i
mean, I am simulating particles in a box, and when a particle goes off
one edge of the box, it comes back on the other side... Easy stuff -
very usesful in molecular dynamics.

<see example code below>

As dist(Maths::dvector& va, Maths::dvector& vb) is an inline function,
it has overall good performance. However if i want to now have a
'generalised' simulation class, my OOP experience tells me that I should
have a 'Space' class that manages calls like 'dist' for any arbritary
derived space system. Thus (via virtual calls) when I call getDistance()
in 'InfiniteSpace' I just call dist(), but when I call it in
'PeriodicRectangular' it takes the 'periodic cube' into account.

My problem is therefore one of performance. When I have a pointer to an
instance of a class that derives from 'Space', the compiler can't inline
the code when getDistance() is called, so performance sux. This goes
away if you make specialised classes, but that defeats the whole point
of inheritance and being generalised.

So my question is, is there a way of giving my simulation a Template
Meta Argument e.g. Simulation<Space s> to avoid the need to have
'virtual' all over the place. I have scratched my head and cant really
think of the best way to do this...

Any Ideas Guys??

Many Thanks For Any Comments,
Jon Rea

-------------
Example Code:
-------------

class Space {
public:
Space(){};
virtual ~Space(){};

inline virtual double getDistance(
const Maths::dvector &va,
const Maths::dvector &vb
) const = 0;
}

class InfiniteSpace: public Space {
public:
InfiniteSpace(){};
virtual ~InfiniteSpace(){};

inline virtual double getDistance(
const Maths::dvector &va,
const Maths::dvector &vb
) const
{
return dist(va,vb);
}
}


class PeriodicRectangular: public Space {
public:
PeriodicRectangular(double allsize):
boxSize(allsize,allsize,allsize),
halfBoxSize(boxSize);

virtual ~PeriodicRectangular(){};

inline void getVector(
const Maths::dvector &va,
const Maths::dvector &vb,
Maths::dvector &vavb
) const
{
vavb.diff(vb,va);

if(vavb.innerdot() < sqrSaveDistance) return;

// find closest image of vb.x;
if(vavb.x > halfBoxSize.x) vavb.x -= boxSize.x;
if(vavb.x < -halfBoxSize.x) vavb.x += boxSize.x;

// find closest image of vb.y;
if(vavb.y > halfBoxSize.y) vavb.y -= boxSize.y;
if(vavb.y < -halfBoxSize.y) vavb.y += boxSize.y;

// find closest image of vb.z;
if(vavb.z > halfBoxSize.z) vavb.z -= boxSize.z;
if(vavb.z < -halfBoxSize.z) vavb.z += boxSize.z;
}


inline virtual double getDistance(
const Maths::dvector &va,
const Maths::dvector &vb
{
Maths::dvector vavb;
getVector(va,vb,vavb);
return vavb.mag();
}
)

Simulate( Space* mySpace )
{
// Do some simulating ...
for( int i = 0; i < xxx; i++ )
mySpace->getDistance( blaA, blaB );
}

 

[Victor Bazarov]

I work with molecular simulations. [..]

I am currently looking at 'Periodic boundary simulations'. By this i mean,
I am simulating particles in a box, and when a particle goes off one edge
of the box, it comes back on the other side... [..]

"It" comes back? Are you sure? I thought the whole point is that when
a particle leaves, it leaves forever, and *another one* comes into the
box to play.

As dist(Maths::dvector& va, Maths::dvector& vb) is an inline function,
it has overall good performance. However if i want to now have a
'generalised' simulation class, my OOP experience tells me that I
should have a 'Space' class that manages calls like 'dist' for any
arbritary derived space system. Thus (via virtual calls) when I call
getDistance() in 'InfiniteSpace' I just call dist(), but when I call
it in 'PeriodicRectangular' it takes the 'periodic cube' into account.

My problem is therefore one of performance. [..]

I don't have an immediate solution for you, especially considering that
I don't have a slightest idea about molecular simulations. But if the
molecule goes outside the box, shouldn't you stop using it in your
simulation and instead use the other one (not "it"), which comes into
the box from the opposite side?

The idea here is that instead of trying to simulate the infinite space
and the infinite number of particles, you take a sample of the space,
limiting it to your box, and only use those particles that are inside
your box at any time.

If you continue trying to track those particles, you're essentially
trying to solve the same infinite space problem using the same number
of particles, only now you make it more difficult to track them.

Of course, it all could be pointless blabber, since I don't really
know squat about molecular simulations.

V

 

[=?ISO-8859-1?Q?Szabolcs_Horv=E1t?=]

So my question is, is there a way of giving my simulation a Template
Meta Argument e.g. Simulation<Space s> to avoid the need to have
'virtual' all over the place. I have scratched my head and cant really
think of the best way to do this...

I'm not sure if this will help ... but take a look at this paper
(especially sections 1.3 and 1.3.3):
http://osl.iu.edu/~tveldhui/papers/techniques/

This is from the same webpage that you mentioned in your last post.

-- Szabolcs

 

[F.J.K.]

So my question is, is there a way of giving my simulation a Template
Meta Argument e.g. Simulation<Space s> to avoid the need to have
'virtual' all over the place. I have scratched my head and cant really
think of the best way to do this...

Any Ideas Guys??

Hi Jon,

Why do you want to use polymorphism? Why not just use a common
functional interface? I'm a bit puzzled, because what you describe
seems to be standard stuff for the most simple usage of templates
conceivable.

As I don't know your surrounding program, it is impossible to tell
whether you want to store the space group in some global variable. In
that case boost::any would be your candidate.

Anyway, have a look at the following which should both be fast and
generic:

#include <iostream>
struct InfiniteSpace {
inline double getDistance(){std::cout << "using infinite space\n";}
};
struct PeriodicRectangular {
inline double getDistance(){std::cout << "using periodic
boundaries\n";}
};
template <typename T>
void Simulate(T& t)
{
std::cout << "Calling function "<<__FUNCTION__<<"\n";
t.getDistance();
};
int main() {
bool use_periodic=true;
InfiniteSpace i;
PeriodicRectangular p;
if (use_periodic)
Simulate (p);
else
Simulate (i);
}

 

[Jon Rea]

Hi Jon,

Why do you want to use polymorphism? Why not just use a common
functional interface? I'm a bit puzzled, because what you describe
seems to be standard stuff for the most simple usage of templates
conceivable.

As I don't know your surrounding program, it is impossible to tell
whether you want to store the space group in some global variable. In
that case boost::any would be your candidate.

Anyway, have a look at the following which should both be fast and
generic:

#include <iostream>
struct InfiniteSpace {
inline double getDistance(){std::cout << "using infinite space\n";}
};
struct PeriodicRectangular {
inline double getDistance(){std::cout << "using periodic
boundaries\n";}
};
template <typename T>
void Simulate(T& t)
{
std::cout << "Calling function "<<__FUNCTION__<<"\n";
t.getDistance();
};
int main() {
bool use_periodic=true;
InfiniteSpace i;
PeriodicRectangular p;
if (use_periodic)
Simulate (p);
else
Simulate (i);
}

The issue is this:

Avove is not quite right for my system because there's a single place
where the particular type of Space is stored: so more like this:

struct Space{
inline double getDistance(){std::cout << "base space\n";}
};

struct InfiniteSpaceuiblic Space {
inline double getDistance(){std::cout << "using infinite space\n";}
};

struct PeriodicRectangular: public Space {
inline double getDistance(){std::cout << "using periodic
boundaries\n";}
};

template <typename T>
void Simulate(T& t)
{
std::cout << "Calling function "<<__FUNCTION__<<"\n";
t.getDistance();
};

int main() {
bool use_periodic=true;
InfiniteSpace i;
PeriodicRectangular p;
Space = the_chosen_space;
if (use_periodic)
the_chosen_space = &p;
else
the_chosen_space = &i;

Simulate( *t );
}

now this wont work because *t is of type Space (no matter what it
actually points to) and so you'll always get the <Space> version of
the template function not the appropriate one (i think) - one'd have
to do a dynamic up cast ..

Am I making sense?

Cheers,
Jon

 

[F.J.K.]

Avove is not quite right for my system because there's a single place
where the particular type of Space is stored: so more like this:

Actually that was what I was referring to in above sentence ;-) At this
point you have two choices. You can use boost::any and its any.type() +
anycast <>() to do the the proper typecasting. I think (but I'm only a
newbie and welcome being corrected) that RTTI and typeid() won't work,
if you don't use virtual functions, which is the whole point.

Or you can keep track of what's really in your variable yourself with
an extra variable like use_periodic and do the typecasting yourself.

Most quantum mechanical or molecular mechanics programs I know of keep
track themselves. While theoretically this allows for runtime errors,
practically it gives you total control over when and where you use
polymorphism (manually, in your if clauses) and when you can't afford
the cost. Even very old QM Fortran programs like Gaussian use
polymorphism in that sense.

now this wont work because *t is of type Space (no matter what it
actually points to) and so you'll always get the <Space> version of
the template function not the appropriate one (i think) - one'd have
to do a dynamic up cast ..
Am I making sense?

Yup. That's exactly what dynamic_casts<> are good for. However in
certain situations, union + static_cast would be equally good, maybe
even better. Ask yourself if there really is a common implementation
base for your Space classes. In the example you gave, there isn't. The
only thing they share is a functional interface, the availability of
the getDistance() function. In that case: don't use inheritance.

Take a look at vector<> and list<>. They both implement functions like
push_back() or size(). Yet there's practically no overlap in
implementation. Consequently the STL isn't so stupid as to define a
class container<> so that vector<T>:container<T> and
list<T>:container<T>. Learn from the mistakes of the early class
libraries which did exactly this and were slow as hell and do it the
STL way

#include <iostream>
struct InfiniteSpace{
double getDistance(){std::cout << "using infinite space\n";}
};
struct PeriodicRectangular{
double getDistance(){std::cout << "using periodic boundaries\n";}
};
union Space {
InfiniteSpace infinite_space;
PeriodicRectangular periodic_rectangular;
};
template <typename T>
void Simulate(T& t)
{
std::cout << "Calling function "<<__FUNCTION__<<"\n";
t.getDistance();
};
int main() {
bool use_periodic=false;
Space space;
if (use_periodic)
Simulate(space.periodic_rectangular);
else
Simulate(space.infinite_space);
}