Template Titbit - A Different Perspective

Oliver Schoenborn's article, "Tiny Template Tidbit", in Overload 47 illustrates some of the wonderful things you can do with templates. I particularly liked the way it described the thought processes Oliver went through when designing some curve-fitting software and his clear explanation of the problem he set himself. Whilst reading an earlier draft of "Tidbit" a slightly different solution occurred to me, which I shared with Oliver before publication. With the deadline approaching there wasn't much time for discussion, so we agreed that I should write up those ideas that seemed interesting, but that were not strictly relevant to Oliver's article.

The "Tidbit" article used a fitCurve() function template to illustrate the techniques of interest:

template <class CoordContainer>
Curve fitCurve(
          const CoordContainer& coords) {
  // ...
  // set coord1 to first Coord in coords
  // set coord2 to second Coord in coords
  Coord c = coord1 + coord2;
  // ...
}

The challenge was to make this template work for containers whose elements are: points, objects that contain a point, or pointers to any of those types.

Oliver called his point class "Coord". The key to his solution was a getCoord() function template that returned the Coord corresponding to a given container element. Different specialisations of getCoord() were defined for different types of element.

// General function template
template <class PntClass> inline
const Coord& getCoord(const PntClass& p)
{ return p.coords; }

// Partial specialization for pointers
// to things
template <class PntClass>
const Coord& getCoord(const PntClass* p)
{ return p->coords; }

// Complete specialization for Coord
inline
const Coord& getCoord(const Coord& p)
{ return p; }

// Complete specialization for pointer
// to Coord
inline
const Coord& getCoord(const Coord* p)
{ return *p; }

This version of the getCoords() function template covers element types that are: Coord, pointer to Coord, something containing a 'coords' data member, and pointer to something containing a 'coords' data member. It doesn't cover elements whose Coord is accessed via a member function. For this, a getData() function template was introduced and the appropriate getCoord() specialisations were re-written in terms of getData().

template <class PntClass>
const Coord& getCoord(const PntClass& v) {
  return getData(v, &PntClass::coords);
}

// partial specialization for pointers
template <class PntClass>
const Coord& getCoord(const PntClass* v) {
  return getData(*v, &PntClass::coords );
}

template <class PntClass>
inline
const Coord& getData(const PntClass& pp,
                     Coord PntClass::* dd) {
  return pp.*dd;
}

template <class PntClass>
inline
const Coord& getData(const PntClass& pp,
                     const Coord& (PntClass::* mm)() const) {
  return (pp.*mm)();
}

This is a neat and, as far as I know, original idea that makes it possible to use either data members or member functions interchangeably in generic functions.

I made three suggestions:

Oliver liked the traits idea, but found it didn't give him what he wanted. And, from his perspective, the other two suggestions just weren't relevant to the problem at hand. After some thought I realised that it was all a question of your point-of-view. As Kevlin Henney often points out, the solution to a problem depends on the context. Oliver had chosen the context of a specific project; I was thinking of a general-purpose library.

Let's look at these ideas in a bit more detail.

Passing a pair of iterators to the fitCurve() function is more general than passing a container, but it's often less convenient for the calling code. An iterator range allows containers without begin() and end() functions to be used (such as arrays) and makes it easy to fit a curve to a subset of the points in a given container. On the other hand, it means passing two parameters instead of one.

For a general-purpose library flexibility is very important. Oliver's project doesn't need that much flexibility, so ease-of-use dominates the trade-off. Context matters here.

The "Tidbit" article shows how to write some generic components that support containers in which the element type is in one of four fairly general categories. In that article Oliver makes the point that this may not be general enough for a library and I agree. What I had in mind was a traits class template that defines a mechanism for accessing the co-ordinates of a point class. The following template and partial specialisation illustrates the idea for points in the X,Y plane.

// The general version of the point
// traits class.
template<typename point_type>
struct point_traits {
  static double x(const point_type& p) {
    return p.x();
  }
  static double y(const point_type& p) {
    return p.y();
  }
};

// A partial specialisation of the point
// traits class for pointers.
template<typename point_type>
struct point_traits<point_type*> {
  static double x(point_type* p) {
    return p->x();
  }
  static double y(point_type* p) {
    return p->y();
  }
};

The general version of the template provides co-ordinate access functions for point classes having member functions x() and y(). The partial specialisation does the same for pointers to classes with x() and y() member functions. For any other type of point class the user would have to define a further specialisation of the traits template. For example:

// A complete specialisation of the point
// traits class for POD_Point.
struct POD_Point { double x, y; };

template<>
struct point_traits<POD_Point> {
  static double x(const POD_Point& p) {
    return p.x;
  }
  static double y(const POD_Point& p) {
    return p.y;
  }
};

The point_traits<> member functions are similar to Oliver's getCoord() functions. They are template functions and their purpose is to extract information from some sort of point object. They are used in the curve fitting algorithm in a similar way. My implementation of a least-squares line algorithm, for example, contains the following lines:

template<typename point_type>
void add(point_type point) {
  double x = point_traits<point_type>::x(point);
  double y = point_traits<point_type>::y(point);
  // ...
}

In principle, the traits member functions could be implemented in terms of another function template, like Oliver's getData(), which would make it possible to handle both data members and member functions out of the box. I prefer not to do this, though. The traits mechanism is general enough to be used with arbitrary types of point object and the inconvenience of having to define a traits specialisation for point classes without the x() and y() access functions is small. In fact, the lack of support for data members could even be seen as an advantage - it should encourage the use of fully-fledged, properly encapsulated point classes. More importantly, though, I wouldn't want the name of a member of a class in the client code to appear in my library functions the way 'coords' appears in Oliver's getCoord() functions. If we were to do this in a general purpose library what would we call the member? 'coords', 'point', 'm_2Dpoint', 'xyPoint_', … The list is endless!

Of course, if you are working on a project that already has to deal with a variety of point classes, some of which use data members, the need for extra traits specialisations could be a nuisance. But, then again, how many different point classes do you want to use in the same application? Not many, I think.

Once again, context matters.

The only curve fitting algorithm I am familiar with is the least-squares line and its 3-dimensional cousin, the leastsquares plane. In fact, the least-squares line can be defined as follows:

A straightforward implementation of this formula in code involves four totals (∑ x_i, ∑ y_i, ∑ x_ix_i and ∑ x_iy_i) and a point count. Old fashioned C++ code might look like this…

// Initial data.
struct Point {double x, y;};
Point point[20] = { ... };
unsigned n = 0;
double x_sum = 0, y_sum = 0,
       xx_sum = 0, xy_sum = 0;

// Accumulate sums for each (x,y)
// point...
for (int i = 0; i < 20; ++i) {
  ++n;
  x_sum += point[i].x;
  y_sum += point[i].y;
  xx_sum += point[i].x * point[i].x;
  xy_sum += point[i].x * point[i].y;
}

// Calculate slope and offset.
double m = (n * xy_sum - x_sum * y_sum)
         / (n * xx_sum - x_sum * x_sum);
double c = (y_sum - m * x_sum) / n;

The algorithm reminded me of the standard accumulate function, so my first thought for a least-squares function in the modern style was modelled on std::accumulate().

// A line in the X,Y plane.
class line;

namespace least_squares {
  // The current state of the least-
  // squares line algorithm.
  struct state {
    state() : n(0), x_sum(0), y_sum(0),
              xx_sum(0), xy_sum(0) {}
    operator line() { ... }
    // ...
    unsigned n;
    double x_sum, y_sum,
           xx_sum, xy_sum;
  };

  // Fit a least-squares line to the
  // (x,y) points in [first,last)
  template<typename iterator_type>
  state fit(iterator_type first,
            iterator_type last,
            state = state());
}

The curve fitting function is given an initial state and returns a new state after accumulating points in the iterator range [first,last). This is what I have called a "function with external state". The default initial state corresponds to no points accumulated.

The state class has a conversion operator that provides an implicit conversion to a line. The combination of external state and conversion operator makes it possible to fit a least-squares line in several stages.

using least_squares::state;
using least_squares::fit;

state algorithm_state =
             fit(point, point + 10);

line best_fit_line =
             fit(point + 10, point + 20,
                 algorithm_state);

Unfortunately, there are several unpleasant features in this design. The state class has public data and an implicit conversion operator. Passing the state into and out of the fit() function by value leads to unnecessary copying. And the simple case of fitting a least-squares line to all the points in a container carries the excess intellectual baggage of the state object.

In thinking about these points I finally settled on a design offering two interfaces that differ in generality and convenience. The less general/more convenient interface is provided as a single template function, least_squares::fit().The more general/less convenient interface is provided as two classes: least_squares::state and least_squares::add_point. The point traits idea has been retained to allow the underlying implementation to adapt to different point classes.

// Summary of least-squares components
// for the X,Y plane.
class point;
class line;
template<typename T> struct point_traits;

namespace least_squares {
  // General interface.
  class add_point;
  class state;
  // Convenient interface.
  template<typename iterator_type>
  line fit(iterator_type first,
           iterator_type last);
}

The state class is at the heart of this new design. It now properly encapsulates its data and provides a minimal interface. It has just two member functions: one to add a point and one to create a line.

// The state of the least-squares
// line algorithm.
class least_squares::state {
public:
  state();
  template<typename point_type>
  void add(point_type point);
  ::line line();
private:
  unsigned n;
  double x_sum, y_sum, xx_sum, xy_sum;
};

The add() function uses the point traits template to extract individual co-ordinates from each point, calculates new values for the four totals and increments the point count.

template<typename point_type>
void least_squares::state::add(
                            point_type point) {
  ++n;
  double x =
          point_traits<point_type>::x(point);
  double y =
          point_traits<point_type>::y(point);
  x_sum += x;
  y_sum += y;
  xx_sum += x * x;
  xy_sum += x * y;
}

The line() function calculates the slope and offset of the leastsquares line and returns a line object. At least two points must have been accumulated before this function is called. The pre-condition is checked using assert(), here; throwing an exception or the use of a policy class might be more appropriate for a real library.

::line least_squares::state::line() {
  assert(n > 1);
  double m = (n * xy_sum - x_sum * y_sum)
           / (n * xx_sum - x_sum * x_sum);
  double c = (y_sum - m * x_sum) / n;
  return ::line(m, c);
}

The state class actually provides everything we need to find the least-squares line for a set of points. On its own, however, it is slightly inconvenient to use. The programmer using the curve fitting library needs to iterate through the points, either by writing a loop or by passing a suitable function object to std::for_each(). The add_point class provided by the library is just the right sort of function object for std::for_each(). Its function call operator simply calls the state:add() function.

class least_squares::add_point {
public:
  add_point(state& s) : algorithm_state(s) {}
  template<typename point_type>
  void operator() (point_type point) {
    algorithm_state.add(point);
  }
private:
  state& algorithm_state;
};

With this convenience class, even the general interface to the least-squares algorithm becomes quite easy to use, as the following example shows.

using std::foreach;
using least_squares::state;
using least_squares::add_point;
// Find the line that best fits
// point[0]..point[19].
state best_fit_data;
for_each(point, point + 20,
         add_point(best_fit_data));
line best_fit_line = best_fit_data.line();

The alternative interface just wraps up these three lines of code in a template function.

template<typename iterator_type>
line least_squares::fit(iterator_type first,
                        iterator_type last) {
  state algorithm_state;
  std::for_each(first, last,
                add_point(algorithm_state));
  return algorithm_state.line();
}

The previous curve fitting example can now be re-written as a one-liner, which I leave as an exercise for the interested reader.

So, my answer to the external/internal state question is, "Both". The general interface uses external state in the form of an add_point function object; the convenient interface keeps the algorithm state entirely within the fit() function. Once again, context matters, and this time the library can not guess which is the more important.

The curve-fitting code presented here illustrates how Object- Oriented Programming and Generic Programming can work together. The least_squares::state class obeys the principles of OOP; the fit() function is a generic algorithm. Together they provide a safe, flexible and efficient software component.

The difference between Oliver's code and mine comes down to the context of the problem. It's like spelling. The dictionary in my version of Microsoft Word says "tidbit" is an error and offers "titbit" as the correct spelling. My Chambers dictionary has both spellings. Context matters!